Understanding the Functional Properties of Lipid Heterogeneity in Pulmonary Surfactant Monolayers at the Atomistic Level

1. Introduction

The integrity of the alveolar gas-blood barrier is crucial for effective gas exchange and health, filtering of undesirable components, and response to inhaled hazard. At the same time, it develops tolerance mechanisms to attenuate immunopathology. The alveoli are continuously exposed to inhaled micro- and nanosized pathogens, which are normally rapidly eliminated with the help of the immune system. Immune responses in the alveoli must be tightly regulated to prevent excessive inflammation and tissue damage. Inappropriate or excessive immune responses cause the development of systemic airway inflammation, as in the Acute Respiratory Distress Syndrome (ARDS) (Ware and Matthay, 2000). ARDS is the major cause of respiratory failure affecting millions of people annually, and it is also a main cause of death in many viral infections, such as in severe acute respiratory syndrome coronavirus 1 (SARS-CoV-1), and in the current SARS-CoV-2 which causes the coronavirus disease 2019 (COVID-19) (Xu et al., 2020).

The alveolar epithelium is coated by pulmonary surfactant, that is a delicate membranous proteolipid film that maintains normal lung function. Pulmonary surfactant forms a monolayer, which lines the alveolar epithelium and is synthesized and secreted by epithelial alveolar type 2 cells. Pulmonary surfactant is the most permeable interface of the human body exposed to the environment, presenting the first respiratory barrier against inhaled foreign matter and microorganisms. It is a lipoprotein complex comprising approximately 90 weight-% lipid of which phosphatidylcholine (PC) is the principal component. Importantly, the pulmonary surfactant is exceptionally rich in dipalmitoylphosphatidylcholine (DPPC) with a high main transition temperature (T_m). Other major lipid components involve PCs with an unsaturated chain, phosphatidylglycerol (PG), and cholesterol (Goerke, 1998). As to protein content, the pulmonary surfactant contains four specific surfactant proteins (SP-A, SP-B, SP-C, and SP-D) of which SP-A and SP-D are innate immune defense proteins, whereas SP-B and SP-C together with the phospholipids are crucial in sustaining the very low surface tension needed to avoid alveolar collapse, oedema, and lack of oxygenation (Possmayer et al., 2010).

Basic cellular and molecular biology research has suggested that early pulmonary surfactant dysfunction contributes to the high morbidity of coronaviruses (Xu et al., 2020). The ability to reduce the surface tension of pulmonary surfactant can be compromised due to decreased concentration of surfactant phospholipids and proteins, altered phospholipid composition, proteolysis, and/or protein inhibition, as well as oxidative inactivation of lipids and proteins (Postle et al., 1999; Taeusch et al., 2005; Zasadzinski et al., 2005; Dushianthan et al., 2014). Variations of these dysfunctional mechanisms have been reported in child and adult patients with ARDS. For instance, pulmonary surfactant metabolism studies in adult ARDS patients showed altered surfactant lipid composition: DPPC content was decreased, whereas the fractions of the surface tension-inactive unsaturated species were increased (Amigoni et al., 2017). Moreover, both the total amounts of PC and PG were decreased (Amigoni et al., 2017). Therefore, the possibility that early surfactant replacement therapy could be beneficial in preventing progression of disease severity is encouraging, given the established harmless profile of the pulmonary surfactant. Even though it is mainly indicated for prematurely born babies, surfactant replacement therapy has also been used in adult ARDS studies (Walmrath et al., 1996, 2002; Gunther et al., 2002; Moller et al., 2003). This method (Echaide et al., 2017), where exogenous surfactant is supplied into the lungs, is currently being tested in clinical trials to treat COVID-19 infected patients that require ventilator support (Grocott, 2020; Lewis, 2020). Accordingly, for the development of efficient surfactant replacements, we require a better understanding of the roles of the pulmonary surfactant components in lung mechanics.

Pulmonary surfactant forms a network of complex biological self-assembling morphologies lining the alveoli. The distinctive structure formed by the pulmonary surfactant is a monolayer at the gas-alveolar epithelium liquid interface. In addition to lateral packing in the monolayer at the liquid-air interface, a fraction of the pulmonary surfactant is also likely folded from the interface into lipid bilayers or multilayers in the aqueous subphase (Baoukina et al., 2007a), acting as lipid reservoirs (Pérez-Gil, 2008; Parra and Pérez-Gil, 2015). This monolayer is repeatedly compressed and expanded during breathing cycles without hysteresis—a property that can only be withheld by a material with very peculiar viscoelastic properties (Andreassen et al., 2010; Steimle et al., 2011). These reservoirs have also been suggested to participate in oxygen transfer from the inhaled air to circulation (Olmeda et al., 2010). It has been shown that the pulmonary surfactant monolayers and membranes exhibit phase behavior that is believed to play roles in lung mechanics; tight packing of lipids with saturated lipid chains promotes its ability to lower surface tension, whereas lipids with unsaturated chains increase pulmonary surfactant fluidity and thus allow for its rapid adsorption to the air–water interface (Bernardino de la Serna et al., 2004, 2013a,b; Casals and Cañadas, 2012). Moreover, the phase behavior is also considered to regulate the functions of the surfactant proteins (Pérez-Gil, 2008; Bernardino de la Serna et al., 2009). However, the biophysical implications of the phase behavior are poorly understood, and a molecular view of the molecular organization would be extremely helpful in understanding the roles of different lipids and surfactant proteins in lung functionality enabled by pulmonary surfactant.

When compressed under the T_m of the corresponding membrane, single-component monolayers transition from a two-dimensional gas phase into a loosely packed and dynamic liquid-expanded (L_E) phase. Further compression takes the monolayer to a brittle liquid-condensed (L_C) phase, which eventually collapses, excreting matter into the aqueous subphase (Kaganer et al., 1999; Baoukina et al., 2008; Brewer et al., 2010). The L_E-L_C transition takes place through the formation of a coexistence phase characterized by a plateau in a surface pressure-area isotherm, which is manifested as observable domains (Klopfer and Vanderlick, 1996; Mansour and Zografi, 2007; Javanainen et al., 2017a).

The behavior of native pulmonary surfactant monolayers is quite different from such model systems. Still, fluorescence and Brewster angle microscopy reported visible domains in monolayers formed from lipid fractions of the pulmonary surfactant (Discher et al., 1999; Bernardino de la Serna et al., 2013a). These domains were suggested to be highly enriched in DPPC (Discher et al., 1999). Domains were also detected in vesicles created from the pulmonary surfactant, and also in the presence of the surfactant proteins (Bernardino de la Serna et al., 2004, 2013a). However, the surface pressure-area isotherms measured for pulmonary surfactant lipid fractions were not found to contain a plateau of any sort (Discher et al., 1999). Thus, the monolayer did not entirely transition to the L_C phase at any pressure, indicating that no L_E-L_C coexistence was present either. Presumably, the observed domains did not present an equilibrium phase. Indeed, the composition of the pulmonary surfactant seems to be adjusted so that the mixture is barely fluid at body temperature (Bernardino de la Serna et al., 2004; Suri et al., 2012), which can lead to distinct behavior of its lipid components. By doping the pulmonary surfactant lipid fraction with additional DPPC, the T_m can be gradually increased, and a coexistence plateau becomes eventually visible (Discher et al., 1999). These observations suggest that the behavior of the pulmonary surfactant may be characterized by transient heterogeneity arising from critical fluctuations (Nielsen et al., 2000), where the behavior of each lipid type is determined by its T_m value.

These research topics are difficult to study at the molecular level due to the limited resolution of available experimental techniques. Molecular simulations are often helpful in resolving such aspects down to the atomistic level (Enkavi et al., 2019). However, previous simulation studies on multi-component monolayers are relatively scarce (Baoukina et al., 2007a,b, 2008, 2012, 2014; Baoukina and Tieleman, 2011) with only a few studies performed using atomistic models (Javanainen et al., 2010; Olzyńska et al., 2020). While previous simulations have generated a lot of insight to better understand the behavior of surfactant monolayers, they have been haunted by the insufficiently accurate description of the physical behavior of the lipids at the interface, that, in turn, was largely due to insufficient quality of water models in describing phenomena at water-monolayer-air interfacial regions (Baoukina et al., 2007b; Vega and De Miguel, 2007; Lamberg and Ollila, 2015; Javanainen et al., 2017a). Moreover, The description of the driving forces behind membrane heterogeneity in often used coarse-grained Martini model is limited (Davis et al., 2013).

To overcome these limitations and provide a detailed and accurate picture of the molecular-level organization of surfactant monolayers, we used extensive state-of-the-art all-atom molecular dynamics (MD) simulations of model systems, validated with experimental Langmuir trough data. The composition of our model monolayers was chosen to match the composition of the protein-free pulmonary surfactant (Bernardino de la Serna et al., 2013a). These model systems were simulated at various compression states and temperatures, matching the conditions used in experiments. By using our recently developed protocol for performing simulations of lipid monolayers at the air–water interface, Javanainen et al. (2017a) and Paananen et al. (2019) we were able to reach quantitative agreement with experimental surface pressure-area isotherms. Our simulations predict the presence of lateral heterogeneity characterized by domain formation, which we confirmed by atomic force microscopy (AFM) imaging. Upon compression, we found the appearance of transient L_C-like domains. Upon further compression these domains were found to aggregate to form large ordered regions. Thanks to the atomistic detail available in the simulations, we were able to draw conclusions as to the physical and chemical properties of the nanoscale domains. Interestingly, simulations predict that the domains are not substantially enriched by any lipid type, indicating that the peculiar viscoelastic properties of the pulmonary surfactant arise from the collective behavior of the mixture rather than from the features of its individual lipid components. Our findings help understand lung mechanics and are valuable in the development of strategies to tackle lung conditions, while also paving the methodological way for future studies of the pulmonary surfactant.

2. Methods

2.1. Atomistic Molecular Dynamics Simulations

Until recently, the poor description of surface tension by the commonly used water models (Chen and Smith, 2007; Vega and De Miguel, 2007) has prevented quantitative simulation studies of lipid monolayers at the air–water interface. However, very recently we demonstrated that a combination of the four-point OPC water model (Izadi et al., 2014) and the CHARMM36 lipid model (Klauda et al., 2010) reproduces experimental surface pressure-area isotherms in single-component 1-palmitoyl-2-oleoylphosphatidylcholine (POPC) and DPPC monolayers (Javanainen et al., 2017a), and the agreement is largely quantitative. In this study, we extend this approach to quaternary lipid monolayers whose composition was chosen to match the composition of the protein-free pulmonary surfactant (Goerke, 1998) as accurately as possible. In practice, our systems contained 60 mol% DPPC, 20 mol% POPC, 10 mol% 1-palmitoyl-2-oleoylphosphatidylglycerol (POPG), and 10 mol% cholesterol. The systems were set up as follows. Two monolayers, separated by a water slab and each containing 169 lipids were first set up at an area per lipid (APL) equal to 50 Ų. Next, these monolayers were either expanded or compressed during a 10 ns simulation to an average APL value of 100 or 40 Ų, respectively, using the MOVINGRESTRAINT and CELL keywords in the PLUMED 2.2 package (Tribello et al., 2014). Structures at a total of 17 APL values were extracted and used as initial structures for the production monolayer simulations in the NVT ensemble for 1 μs each. All simulations were performed at both 298 and 310 K, and additional repeats as well as larger simulations of up to 3,042 lipids were performed to provide further validation for the results.

The CHARMM36 model for phospholipids (Klauda et al., 2010) and cholesterol (Lim et al., 2011) was used together with the four-point OPC water model (Izadi et al., 2014). The simulations were performed with the version 5.1.x of the GROMACS simulation package (Abraham et al., 2015), and the recommended simulation parameters for the CHARMM36 force field (Lee et al., 2015) were used to reproduce realistic monolayer behavior (Javanainen et al., 2017a). Details on the simulation setups and on the simulation methodology are provided in the Supplementary Material.

2.2. Analyses of MD Simulations

2.2.1. Surface Pressure-Area Isotherms

Monolayer behavior was characterized and compared to experiments using surface pressure-area isotherms. Monolayer surface pressure π at an APL of A was calculated as π(A) = γ₀ − γ(A), where γ and γ₀ are the surface tensions of the monolayer-covered and plain air-water interfaces, respectively. Here, the surface tensions are obtained from the pressure components along the monolayer plane (P_L = (P_xx + P_yy)/2) and normal to it (P_N = P_ZZ) as γ = L_z × (P_N − P_L)/2, where L_z is the simulation box size normal to the monolayer plane, and a factor of two indicates the presence of two air–water interfaces in the simulation box. The values were obtained with the gmx energy tool provided with GROMACS (Abraham et al., 2015).

2.2.2. Monolayer Isothermal Compressibility

The values for the monolayer isothermal compressibility modulus ($C_s^{-1}$) were calculated from the surface pressure-area isotherms for both the simulated and measured systems as $C_s{}^{-1}=-A\frac{d\pi }{dA}$, where A is the area per molecule at a given pressure and $\frac{d\pi }{dA}$ is the slope of the surface pressure-area isotherm. The derivatives were calculated using Matlab from a piecewise smoothing spline fitted to the isotherms. Smoothing parameters of 0.01 and 0.25 were used for the fit to the simulated and measured data, respectively. With these values, the fit was well within the error bars of the original data, while allowing to avoid unphysical jumps in the compressibility modulus.

2.2.3. Detection of Monolayer Domains

Domains packed like in the L_C phase were detected by clustering lipid chains and cholesterols based on their packing in the monolayer plane. The 10th carbons in the lipid chains and the C14 atom in cholesterol ring were included in the clustering that used the DBSCAN algorithm (Ester et al., 1996). The chosen atoms in the lipid chains capture well the hexagonal packing in ordered structures (Javanainen et al., 2017b), and the C14 carbon of cholesterol resides at the same depth. This carbon is part of both the five-member (D) and six-member (C) rings. For DBSCAN, we used a cut-off of 0.71 nm and a minimum neighbor count of 6. The cut-off was set to the distance at which the first minimum appears in the radial distribution function of the clustered atoms. The clusters were considered to be part of the L_C-like domain. This clustering was performed on conformations separated by 1 ns. The L_C-like fraction, number of individual clusters, and the largest cluster size were extracted for each conformation. All these quantities were then averaged over the trajectory independently for both monolayers in the simulation system. Their average values are reported, whereas the differences in the values extracted for the two monolayers serve as error estimates. Pattern-matching was used to find all residence events in the L_C clusters. These times were histogrammed and fitted with a power law with an exponent of b.

2.2.4. Contact Fraction

Lateral demixing of lipids with unsaturated and saturated chains was quantified by the contact fraction. Following Domański et al. (2012), we defined the contact fraction as f_mix = c_US−S/(c_US−US + c_US−S), where c_US−S is the number of contacts between the lipids with an unsaturated chain (POPC and POPG) and the lipids with saturated chains (DPPC). For a contact, the phosphorus atoms of the lipids had to be within 1.1 nm from each other. We extracted contact data using the gmx mindist tool every 1 ns and calculated the average values over the last 500 ns of the simulations for both monolayers. We plot the mean of these two average values, whereas their difference serves as an error estimate.

2.2.5. Cholesterol Clusters

Cholesterol molecules were considered to be part of a cholesterol cluster if they were in contact with at least one other cholesterol molecule. We used a cut-off of 0.94 nm based on the first minimum in the cholesterol-cholesterol radial distribution function. All residence events were used in the calculation of the probability distribution of cholesterol cluster sizes. The distances were measured from the centers of mass of the cholesterol molecules.

2.2.6. Diffusion Coefficients

Diffusion coefficients were calculated to characterize monolayer dynamics and to detect changes in monolayer packing. The diffusion coefficients were extracted from center-of-mass (COM) trajectories. The motion of lipids with respect to the movement of the monolayer as a whole was analyzed to eliminate possible artifacts due to monolayer drift. The diffusion coefficients were extracted from linear fits to time- and ensemble-averaged mean-squared displacement at lag times between 10 and 100 ns. Two values were extracted—one from each monolayer—and averaged, and their difference served as an error estimate. The GROMACS tool gmx msd was used, and the COM trajectories were generated using gmx traj.

2.2.7. Monolayer Thickness

Monolayer thickness was used to couple AFM height profiles of heterogeneous membranes to the lateral packing within the domains. This thickness was estimated from density profiles calculated for each lipid type along the monolayer normal (z) using the GROMACS tool gmx density. These profiles were aligned at the phosphorus peak of DPPC, and the lipid was set to begin and end at z values where its density crossed 5% of its maximum value.

2.2.8. Lipid Chain Tilt

Lipid chain tilt was used to characterize persistent L_C-like packing in the monolayers. The tilt angle of lipid chains was calculated as the angle between the z axis and the vector joining the 1st and 16th carbons in the fatty acid chains of phospholipids using the GROMACS tool gmx gangle. The angle distributions were averaged over both chains and both monolayers, and the distributions were fitted with a Gaussian. The location of its maximum and the variance were used as the mean tilt angle and its error estimate, respectively.

2.3. Langmuir Trough Measurements

2.3.1. Langmuir–Wilhelmy Compression Isotherms

By means of a specially designed ribbon Langmuir–Wilhelmy trough (NIMA Technology, UK), compression isotherm assays were performed and surface pressures as a function of molecular area were obtained at constant temperature, as described previously (Dohm et al., 2010). The lipid mixture was identical to that employed in the simulations (see above). The employed Langmuir–Wilhelmy trough has a maximum area of 312 cm² and a minimum of 54 cm² and instead of the canonical rigid barriers uses a continuous Teflon-coated ribbon that by moving symmetrically reduces the available area in the aqueous surface, and therefore compresses the confined sample to carry out the compression isotherm. The pressure is recorded using an electronic pressure sensor and a piece of cellulose and employing the Wilhelmy plate technique, with an estimated error of ±1 mN/m among different isotherms, that were assessed at least in triplicate at 298 K. The compression isotherm displayed in Figure 1 shows the first compression of 44.1 nmol of a DPPC/POPC/POPG/Chol mixture of mimetic surfactant monolayer, spread on a saline solution containing 150 mM NaCl. After placing the lipid mixture onto the aqueous subphase surface, the generated lipid layer was left to equilibrate for 10 min, so the solvents could evaporate, and the lipids organize and disperse at the interface. Thereafter, the narrowest area was achieved by compressing the monolayer at a constant speed of 150 cm²/min. This compression speed was set to reflect the functional quality of the assessed surfactant, by recapitulating as much as possible the dynamic context of native pulmonary surfactant in the human organism, and allow reaching the highest surface pressure possible.

Figure 1. Surface pressure—area isotherms from simulations and experiments, with AFM images and simulation snapshots taken from selected surface pressures at 298 K. Scale bar 1 μm. At 25 mN/m, (3) one can identify small irregularly shaped domains with visible nanostructures inside the domains. At 35 mN/m, (2) the domains start to fuse and lose their rounded shapes. In the collapse plateau at 45 mN/m, (1) the domains are no longer visible in the AFM images, and some lipids are squeezed out from the monolayer film. Clearly visible transient or persistent heterogeneity can also be observed in the simulation snapshots (4–7) of monolayers at different surface pressures. Here, (4) was extracted from a 4-fold larger monolayer than (5–7). Only the data points below the experimental collapse pressure are shown in the isotherm, whereas all simulated points with all of the independent repetitions are shown in Supplementary Figure 1.

2.3.2. Surface Dilatational Rheology

We again used a lipid composition identical to the one in simulations (see above). All these lipids were obtained from Avanti Polar Lipids (Alabaster, AL). Lipids were spread in 2 mM chloroform solution on the surface of a deionized water subphase (Milli-Q, Millipore, Bedford, MA) in a Langmuir trough (KSV Minitrough, Espoo, Finland) with a platinum Wilhelmy plate. The dilatational rheology of the monolayer was studied using the oscillating barrier method. The film was first compressed to the desired surface pressure, after which sinusoidal area compressions with an amplitude of 1% were performed at a frequency of 10, 50, 100, or 200 mHz and the changes in surface pressure in response to the oscillations were recorded. This allowed us to determine the surface dilatational modulus, (E = −dΠ/d ln A), whose real part represents the elastic modulus, and imaginary part the viscous modulus. In the limit of slow oscillations, the elastic modulus corresponds to the compressibility modulus introduced earlier. Each measurement was repeated two times, and an average dilatational modulus was calculated over frequency, since the magnitude of the dilatational modulus was found to be constant with respect to frequency, although the ratio of compressibility to viscous moduli decreased upon increasing frequency. All measurements were conducted at room temperature (297 K).

2.4. Atomic Force Microscopy Imaging

2.4.1. Langmuir–Blodgett Transferred Monolayers

Lipid monolayers were transferred onto mica as described in Brown et al. (2011). The transfer ratio was equal to 1, i.e., no compression or expansion of the monolayer took place during the transfer. For the monolayer preparation, the sample was spread at 298 K onto the air–water interface until the minimum surface pressure of ~0–1 mN/m was observed. After 10 min of monolayer equilibration, the film was compressed until the desired surface pressure was reached (from 15 to 45 mN/m, in 5 mN/m steps), at a compression speed of 50 cm²/min. Before the transfer was started, the film was again equilibrated for 5 min at constant pressure. The monolayers were finally deposited in a freshly cleaved muscovite mica substrate (Plano GmbH, Wetzlar, Germany) that had been previously submerged. The lifting device—to which the mica substrate was fixed—was raised in the vertical plane out of the buffered aqueous subphase at a speed of 10 mm/min at constant pressure. In all experiment modalities three independent experiments were carried out, as a minimum, and up to 10 images were taken and analyzed.

2.4.2. Atomic Force Microscopy

Langmuir–Blodgett supported monolayers' topographical images were taken using an atomic force microscope (JPK NanoWizard, JPK Instruments, Berlin, Germany), employing in both cases Silicon-SPM cantilevers (Nanosensors, NanoWorld AG, Neuchatel, Switzerland). The AC mode in air was selected for monolayers. The scan rate was ~1 Hz for all AFM images. At least three different supported monolayer systems were assessed, and each sample was imaged on a minimum of three different positions. Image processing of AFM data was done using the SPIP software package as in Bernardino de la Serna et al. (2013a) (Image Metrology, Hørsholm, Denmark).

3. Results

3.1. Monolayers Display Lateral Heterogeneity Without a Coexistence Plateau

Surface pressure-area isotherm is a key quantity representing monolayer behavior at the air–water interface, and it is readily extracted from both Langmuir trough measurements and computer simulations. These isotherms for the quaternary pulmonary surfactant lipid monolayers at 298 and 310 K are shown in Figure 1.

The isotherms from experiments and simulations are in nearly quantitative agreement. No plateau indicating L_C/L_E coexistence was visible in the experimental isotherm, and the isothermal compressibility modulus extracted from these isotherms (see Supplementary Figure 16) remained well below the value of 250 mN/m, typically assigned to the L_C phase, suggesting that the monolayers remained in the L_E phase. Still, our AFM measurement displayed heterogeneity in lipid packing. The AFM images taken at 298 K (Figure 1, panels 1–3, and Supplementary Figures 2–5) revealed surface pressure-dependent formation of a heterogeneous monolayer with thinner and thicker regions. Starting from low surface pressure (higher APL), at 25 mN/m (Figure 1, panel 3) the observed domains were rounded yet irregular in shape. The insides of these domains are heterogeneous with visible nanostructures, indicated also by our height profile studies of the domains (see Figure 2, Supplementary Figures 6–9). Similar heterogeneity was also evident in the simulations (Figure 1, panels 4–7), where domains were readily seen at surface pressures above 25 mN/m, or at APL below 65 Ų at 298 K. No L_C/L_E coexistence plateau was observed in the simulated isotherms either (Figure 1, Supplementary Figure 1), and the isothermal compressibility modulus also followed the experimental values reasonably well (Supplementary Figure 16). Moreover, independent replica simulations at selected APLs, as well as additional simulations of larger monolayer systems confirmed that the initial monolayer configuration and the finite-size effects did not affect the extracted values of surface pressure (Supplementary Figure 1) nor the presence of heterogeneities. Still, the lipid organization in the larger systems is most likely affected by their initial configurations, as lipids do not have time to sample the whole monolayer during the simulation timescale. Therefore, we did not consider the larger monolayer systems in further analyses. The movies at doi: 10.6084/m9.figshare.12612317 as well as the additional snapshots of the simulated monolayers at 298 and 310 K given in Supplementary Figures 10, 11, respectively, clearly demonstrate heterogeneity in lateral organization.

Figure 2. The AFM height profile at 35 mN/m at 298 K and the average monolayer thickness calculated from the simulations. (A) The AFM height profile results show a difference of approximately 0.7 nm at the border of the two phases at the selected positions. By considering the average height of the phases and averaging over all boundaries captured in the profile, the average thickness difference is similar, at 0.52 ± 0.15 nm. (B) The path of the measurement. Two points marked with the blue and red triangles indicate selected positions at the border of the L_C- and L_E-like regions. The scale bar at the lower right corner equals to 1 μm. (C) The average thickness of the monolayer in the simulations. The dashed line indicates the position of maximum of phosphorus density of DPPC. The lines indicate the span of lipids from tails to head groups. At each extreme, the line shows the location where the lipid density drops below 5% of its maximum value.

Based on our simulations at 298 K (Figure 2C), the monolayers at APLs equal to 47.5–50 Ų (corresponding to the very packed state with major regions in an L_C-like arrangement) and 90–100 Ų (corresponding to the L_E phase) extend approximately 2.5 nm and 1.7–1.8 nm toward the air phase, respectively, calculated from the plane defined by the DPPC phosphorus atoms (see section 2). This difference of 0.6–0.7 nm matches well the differences in height profile observed in our AFM experiments at surface pressures below monolayer collapse (Figures 2A,B, Supplementary Figures 6, 7), indicating that the measured thicker and thinner regions in the monolayers would correspond to L_C- and L_E-like regions, respectively. The thickness difference of these two regions in simulations is highlighted in Supplementary Figure 17.

The extent of the L_C-like region grows as the surface pressure of the monolayer is increased, as would also be expected during a proper phase coexistence region (Javanainen et al., 2017a), although the surface pressure does not remain constant here based on the isotherms shown in Figure 1. From AFM figures, at 35 mN/m (Figure 1, image 2) the L_C domains seem to fuse together to form more irregular shapes, but the overall heterogeneity is clearly visible. Here, the coalescence of the domains was likely also limited by the slow diffusion within the compressed monolayers. Indeed, our simulations predict that compression from a state with APL equal to 80 Ų to the one with APL equal to 50 Ų slows lipid diffusion by 1–2 orders of magnitude (Figure 3A). The nanostructures within the domains correspond to the L_C-like ones measured at 25 mN/m (see Supplementary Material). At 45 mN/m (Figure 1, image 1 below APL 55 Ų) and 55 mN/m (Supplementary Figure 5) we have already reached the monolayer equilibrium collapse plateau seen in the experimental surface pressure-area isotherms (Mansour and Zografi, 2007). It should be noted that the monolayers in the alveoli are expected to reach surface tension values close to 0 mN/m (corresponding to surface pressure of ~70 mN/m). However, in experiments this is only achieved by very rapid compression (Crane and Hall, 2001), or by the use of bubble surfactometer methods (Schurch et al., 1989). In the collapse region, the domains are completely merged with a continuous network of small holes, forming a sponge-like collapsed structure at the boundary of the air-liquid interface. Further on, the protruding regions seen in the AFM images at surface pressures above the collapse plateau are collapsed material, excluded from the interfacial film into the water phase. This is verified by the height profiles (Supplementary Figure 9), which show height differences of up to 10 nm, which is clearly above the difference observed in the simulations between monolayers at different APLs (Figure 2C). Increasing the temperature to 310 K in the simulations renders the monolayer more fluid and the heterogeneity is visible at surface pressures above 30 mN/m, or APL below 60 Ų (Supplementary Figure 11).

Figure 3. (A) Diffusion coefficients D of the lipids in the monolayers. Top row: The diffusion coefficients, shown in linear scale as a function of APL, show no signs of a phase transition at either studied temperature. Bottom row: D, in logarithmic scale, shown as a function of 1/APL displays two slopes with a crossover at APL≈60 Ų (298 K) or ≈55 Ų (310 K), highlighted with a dashed gray line. (B) Tilt angles of the phopsholipid chains in the monolayers at 298 K. The shaded areas, bordered by dim lines, show the standard deviation. Data for 310 K is provided in Supplementary Figure 12. (C) Dilatational modulus of the monolayer as a function of surface pressure at 297 K.

Concluding, both AFM measurements and computer simulations revealed lateral heterogeneity in a broad range of APL values, despite the lack of a plateau in the surface pressure-area isotherm. Together, simulations and AFM imaging suggest that the structure of these domains resembles the L_C phase.

3.2. Condensed Nanodomains Affect the Structural and Dynamic Properties of the Monolayer

As no coexistence plateau is observed in the isotherms obtained from experiments and simulations, it is worth asking whether the heterogeneity observed by AFM and simulations manifests itself in some other monolayer properties. To this end, we measured the dilatational modulus of the monolayer at various surface pressures using the oscillating barrier approach. As shown in Figure 3C, we observe little change in this modulus up to a surface pressure of 30 mN/m. This signals that the monolayer area can be changed relatively easily, possibly due to a rapid reorganization of lipids in the membrane plane. However, beyond 30 mN/m—corresponding roughly to APL equal to 60 Ų based on the isotherms in Figure 1—the modulus grows rapidly. This indicates that there are structural changes in the monolayer, perhaps due to the formation of continuous or percolating ordered regions.

On the simulation side, we analyzed the diffusion of lipids at different APLs. The diffusion coefficients, D, are shown in the top row of Figure 3A. The curves demonstrate what seems to be a continuous and smooth increase in diffusion coefficient with increasing APL. However, as shown in the bottom row of Figure 3A, a plot of log(D) vs. 1/APL does not reveal a single slope, as expected for monolayers without a plateau in the pressure–area isotherm, such as DLPC at 294 K (Peters and Beck, 1983; Mangiarotti et al., 2014). Instead, we observed two distinct slopes with crossovers at 60 and 55 Ų for DPPC at 298 and 310 K, respectively. The presentation of the diffusion data used herein is based on the idea that the diffusion of lipids depends on their degree of packing. This presentation resembles that used to interpret diffusion in terms of the free volume theory (Almeida et al., 1992). However, due to various issues in applying this theory to diffusion in lipid monolayers (Javanainen et al., 2010), we refrain from interpreting the diffusion data further within this framework. Yet, we acknowledge a clear change in the trend of the APL-dependence of diffusion coefficients.

Let us move on to a structural property extracted from the simulations. The tilt angle distributions at 310 K, shown in Supplementary Figure 12, show a continuous increase in the chain tilt upon increasing area. However, as show in Figure 3B, at 298 K and between APLs of 45 and 60 Ų, the chains remain tilted at an angle of ~25°. A similar tilt angle of ~25° was also reported for DPPC monolayers in the L_C phase using both X-ray diffraction (Lee et al., 2002) and vibrational sum frequency generation spectroscopy (Ma and Allen, 2007). This indicates that at 298 K, the monolayer retains regions of L_C-like packing with a characteristic chain tilt until relatively large APL values, and this structural feature is coupled to a diffusion mode that differs from that observed at larger APL values. It is worth noting that the long-range orientational order of the tilted L_C-like domains (Moy et al., 1986) is also clearly evident in Figure 1.

Concluding, in both experiments and simulations we systematically observe a a sudden change in monolayer properties at a certain compression level, despite the isotherms showing no coexistence plateau.

3.3. Condensed Domains Display Little Enrichment in DPPC and Cholesterol

With the presence of L_C-like domains evident in both simulations and experiments, we next use the simulation data to study the composition of the domains. The leftmost panels in Figure 4A show the fraction of phospholipid chains and cholesterols in L_C-like domains obtained from clustering analyses (see section 2). Examples of the results of the clustering algorithm at different APL values are shown in Supplementary Figure 14. The data are resolved by lipid type and shown in a cumulative manner so that the total L_C-like fraction can be read from Figure 4A. At low APL values, almost the entire monolayer is part of the L_C-like region, yet the excluded lipids can still be readily compressed, resulting in a compressibility modulus that is lower than values associated with the L_C phase. At the other end, at large APL only few transient clusters of a few lipids are observed. Still, more than half of the phospholipid chains and cholesterol molecules are part of an L_C-like domain up to APL values of ~65 Ų. These results naturally depend on the parameters used for the clustering algorithm, as demonstrated in Supplementary Figure 15. Still, our choice of cutoff based on the first minimum in the radial distribution function of the clustered particles (see section 2) provides results at the extreme APL values that are intuitively expected. Interestingly, temperature has little effect on the calculated L_C fractions, since very similar domain compositions are extracted from the simulations at room and body temperatures.

Figure 4. (A) The extent and composition of the L_C-like regions. On the left, the fraction of phospholipid chains and cholesterols in the L_C-like regions is shown for all lipid types. These curves sum to the total fraction of the L_C-like region. The error bars show the differences calculated for the two monolayers in the simulation system. On the right, the fraction of lipids that make up the L_C-like regions at different APLs, i.e., in this plot the curves always sum to 1. The dashed lines show the values expected based on the overall monolayer lipid composition. (B) Mean contact fraction during the last 500 ns between DPPC and the lipids with unsaturated chains (POPC and POPG). Ideal mixing is highlighted by a dashed line at 0.66, whereas values smaller than that indicate a preference to demix. (C) Left: 2-dimensional probability distributions of cholesterol molecules residing in a cluster with at least one other cholesterol at all simulated APLs at 298 K. Note that the color bar is in logarithmic scale. Right: The sum of probabilities of n > 1, i.e., the probability for cholesterol to be in a cluster with size larger than one molecule. Data for simulations at 310 K are given in Supplementary Figure 13.

It is clear from Figure 4A that DPPC with its two saturated chains makes up the majority of the L_C-like regions, regardless of the surface pressure and temperature. This is natural, since DPPC is also the most prevalent lipid type in our simulations. For a more detailed look into the composition of the L_C-like domains, we plot the fraction of lipids in these regions in the right panels of Figure 4A. Here, the dashed lines highlight the overall fractions of the different components in the system. At very low APL, essentially the entire monolayer is detected by the clustering algorithm to be part of a L_C-like region, and therefore the fractions of lipids in this region agree with their fractions in the system. However, at intermediate APLs between 50 and 80 Ų, DPPC is present in the L_C-like domain slightly more than expected by its fraction in the system. The fraction of POPC within the L_C-like regions is fairly constant over all APLs, whereas those of POPG and cholesterol change more noticeably while moving toward larger APLs; the fraction of POPG decreases while that of cholesterol increases in the shrinking and more transient L_C-like regions. Indeed, at large APLs, cholesterol seems to be a key player in the formation of condensed lipid clusters.

It seems somewhat counter-intuitive that the L_C-like domains have so little preference for any lipid type. However, stable L_C and L_E phases exist in monolayers made of a single lipid type (Javanainen et al., 2017a). We analyzed whether the lipid types underwent lateral demixing during the simulation. The mean contact fractions between lipids with unsaturated and saturated chains in all simulations are shown in Figure 4B. The dashed line at 0.66 represents ideal mixing of DPPC with the lipids with unsaturated chains (POPC and POPG). It is evident from Figure 4B that there is generally little preference for lipid demixing, yet at smaller APL values there is a small tendency for DPPC lipids to reside with other DPPC molecules. The APL values below which lipid mixing deviates from ideal behavior are again close to 55 and 60 Ų, where rapid changes in the APL-dependence of many quantities were observed in the last section. Curiously, temperature has again very little effect on this demixing. Finally, it is worth pointing out that we also studied the time evolution of the contact fractions of the systems at different APL values but found no systematic tendency for any of the monolayers to demix during the simulation time.

We also studied whether cholesterol has a tendency to cluster in the monolayers. The probability of a cholesterol to be in a cluster of a given size is shown in the left panel in Figure 4C for the simulations at 298 K, whereas the data for simulations at 310 K are given in Supplementary Figure 13. The right panel shows the sum of probabilities with n ≥ 2, corresponding to the probability that a cholesterol molecule is part of a cluster with at least one other cholesterol molecule.

It is evident from Figure 4C that all distributions peak at 1, indicating that most cholesterol molecules are on average not in contact with other cholesterol molecules. Still, the probability of a cholesterol to be in a cluster with size larger than one is ~20% in all monolayers (right panel in Figure 4C). Curiously, at APLs 65 Ų and larger, clusters of more than 5 cholesterol molecules sometimes appear. Still, the probability of cholesterol clusters larger than one shows no clear dependence on APL, although not surprisingly, the chance of finding cholesterol clusters is highest in the monolayers with the smallest APLs.

To analyze the distribution of the domains in the monolayer, we show in Figure 5A the number of condensed clusters as well as the fraction of the system occupied by the largest cluster. Between APLs of 45 and 60 Ų the size of the largest cluster drops steeply from a coverage of ~80% of the lipids to a mere ~10%. At the same time, the number of L_C clusters increases from a few to 15–20, and these maxima are reached at an APL equal to 60 Ų (298 K) or slightly smaller (310 K). These APL values agree well with those at which changes in trends as a function of APL were observed in many properties in the previous section. The melting of the largest cluster is coupled to the formation of many smaller clusters. Upon an increase in APL beyond 60 Ų, the number of clusters begins to decrease again, until only a few are detected in the pure L_E phase. Moreover, these few are very small, since the largest one occupies only 2% of the lipids, or less, indicating that it is likely a false positive due to density fluctuations and the used clustering algorithm. It is also noteworthy that at 298 K and between 45 and 70 Ų, there are slightly fewer domains as compared to the case at 310 K. However, the largest domain at 298 K is respectively larger in this interval of APLs.

Figure 5. Spatial and temporal extent of the L_C-like regions. (A) The numbers of clusters in monolayers at 298 and 310 K are shown in green and blue, respectively, whereas the largest clusters in these systems are indicated red and orange, respectively. (B) The exponents b of the fits in the form of P(τ) = a × τ^b to the residence time distributions in (C). The error bars show the differences calculated for the two monolayers in the simulation system. (C) The lifetime distributions of the L_C-like domains for both 298 K (top) and 310 K (bottom). The solid lines show the power law fits.

To evaluate whether the observed domains are transient or long-lived, we analyzed the distributions of residence time of the lipid chains in the L_C-like regions. The residence time distributions for selected monolayers are shown in Figure 5C. These distributions clearly show power law scaling, i.e., the probability P of a binding time τ follows P(τ) = a × τ^b. In the log-log scale, we extracted exponents b from a linear fit, and these exponents are shown in Figure 5B. The larger in absolute value of b, the faster the distribution decays, i.e. the less lipids reside for extended times in the L_C-like clusters. The exponents decrease upon increasing APL indicating that the residence times of lipids in the L_C-like regions decrease. At an APL of 50 Ų, there are more than 100 events where a lipid stays in the L_C-like domain for more than 100 ns, and some lipids stay there throughout the entire simulation. Interestingly, the exponents start to rapidly decay at APLs of 60 Ų (298 K) or 55 Ų (310 K), indicating that long-lived domains are replaced by more dynamic heterogeneity. These APL values, again, agree well with those at which significant changes are observed in other monolayer properties (see previous section). At APLs smaller than 55 Ų, the distributions at the two temperatures are very similar due to the large extent of the L_C-like regions, whereas there are clear differences in the APL range between 55 and 75 Ų. Namely, the residence times for monolayers at 298 K are drastically longer. At areas above 75 Ų, the L_C-like clusters are extremely dynamic and the longest residence times do not reach much over 10 ns regardless of the temperature.

4. Discussion

We first extracted surface pressure-area isotherms from simulations and experiments. As these isotherms did not display a coexistence plateau, we performed AFM experiments that demonstrated the presence of sub-micrometer sized domains (Figure 1, Supplementary Figures 2–5). The height profiles (Figure 2B, Supplementary Figures 6–9) and the thickness profiles obtained from our simulations (Figure 2C) indicated that the domains resemble the L_C phase with hexagonal tail organization, whereas the remainder of the monolayer is in the L_E phase. A detailed investigation of the atomistic simulation data on the condensed domains showed that at large APL (>60 Ų) these domains consisted of transient, isolated islands. As long as the condensed domains are “cushioned” by the surrounding L_E phase, the macroscopic properties of the monolayer such as dilatational/compressibility modulus resemble those of the L_E phase. Here, the domain morphology likely plays a minor effect. Interestingly, a similar observation emerged from a recent coarse-grained simulation study on the elastic properties of lipid bilayers (Eid et al., 2020); it was shown indeed that, in the presence of phase separation, bilayer elastic properties (including both the bending modulus and the area compressibility modulus) are close to the values calculated for the softer component. The present investigation extends the previous conclusion to monolayers not showing proper phase separation, but still presenting structural heterogeneity in the form of dynamic nano-sized domains. At smaller APL (<60 Ų) the condensed domains become stable and form a continuous meshwork, which explains the sudden increase in rigidity of the monolayer indicated by the change in the dilatational modulus in Langmuir trough measurements (Figure 3C). The shift from transient to stable condensed domains was also accompanied by a change in the lateral diffusion dynamics in the monolayer, as well as in the lipid mixing properties (Figure 4B) and in the lifetimes of the L_C-like domains (Figure 5).

Stable phase separation has also been reported, employing fluorescence microscopy and AFM, in pulmonary surfactant extract bilayer (Bernardino de la Serna et al., 2004) and monolayer membranes (Bernardino de la Serna et al., 2013a), whose behavior was similar when the monolayer was compressed to a pressure of ~30 mN/m. The monolayer phase behavior is naturally more complex of the two, as it is strongly affected by lateral pressure, as highlighted by AFM experiments (Bernardino de la Serna et al., 2013a). Importantly, the AFM images measured for monolayers from lung surfactant extracts [compare Figure 1 with Bernardino de la Serna et al. (2013a)] are very similar to those obtained for our synthetic monolayers. This suggests that our quaternary lipid mixture captures the central features of the phase behavior of the entire complex pulmonary surfactant, suggesting a minor role of surfactant proteins in this respect. Lateral heterogeneity plays a major role in lung mechanics, as pulmonary surfactant needs to possess certain viscoelastic properties in order to decrease the surface tension to low levels upon compression, and at the same time maintain its ability to spread rapidly at the air–water interface, and fold away from the interface to form lipid reservoirs in the aqueous subphase when the lateral pressure exceeds a certain threshold.

The role of DPPC in the pulmonary surfactant has been pictured to reduce surface tension exclusively by forming condensed domains (Discher et al., 1999). This conclusion was drawn from experiments on calf pulmonary surfactant extracts that displayed large flower-shaped DPPC-rich L_C-like domains. However, these experiments were performed at 293 K, i.e., 17 K lower than the physiologically relevant body temperature, which can have a major effect on phase behavior near a phase transition. In contrast, our simulations indicate that at room and physiological temperatures, DPPC is only slightly enriched in the L_C-like domains. Based on such results, we hypothesize that the role of DPPC is to act—together with cholesterol—as a nucleation center for the formation of domains into which lipids with unsaturated chains can also merge. This is further evidenced by our analyses on domain compositions and lifetimes, which suggests that lipids in monolayers have little tendency to undergo phase separation based on the saturation level of their acyl chains (Figure 4).

In terms of fluidity, POPG has been suggested to be a key player (Hook et al., 1984), in addition to its role in protein-lipid interactions (Pérez-Gil, 2008; Liekkinen et al., 2020). We observed that bulk L_E regions are indeed enriched in POPG. Still, the fluidity is not due to POPG alone, since all the lipid components display fairly similar diffusion coefficients across APL values (Figure 3). Additionally, we observed that cholesterol has a slight tendency to induce ordered clusters at large APLs, whereas the cholesterols themselves show little tendency to cluster together.

The fact that the spatial heterogeneity was dynamic and dictated by different packing (“physical separation”) instead of the nature of the lipid species (“chemical separation”), suggests that the pulmonary surfactant has specific and collective viscoelastic properties that cannot be derived from the behavior of its components in a straightforward manner. These dynamic properties will be further clarified in our future work. This picture is drastically different from that obtained in similar mixtures using the coarse-grained Martini model: Baoukina et al. observed a plateau in the pressure-area isotherm for a monolayer that underwent “chemical” separation into very distinct phases (Baoukina et al., 2012). Thus, the absence of a plateau in experimental isotherms is in line with the domains visible in AFM images not being compositionally very different from the bulk parts of the monolayer. Moreover, our earlier studies (Bernardino de la Serna et al., 2009, 2013a) have indicated similar “physical separation” of lipids into observable co-existing L_E and L_C-like regions that have matching diffusion properties, indicating that the chemical compositions of the phases can remain similar in heterogeneous systems. Our simulations also indicate that the behavior of the monolayers is in general very similar at room and body temperature at the studied surface pressures below 50 mN/m. This is not very surprising, as both temperatures fall below the T_m of DPPC. In addition to these well-balanced viscoelastic properties of the pulmonary surfactant, the observed heterogeneity likely plays a role in regulating the function of surfactant proteins (Pérez-Gil, 2008; Bernardino de la Serna et al., 2009).

Computer simulations have become an indispensable tool in biological soft matter research due to their ability to probe small time and length scales (Enkavi et al., 2019). Traditionally, reproducing surface pressure-area isotherms in classical molecular dynamics simulations has been challenging, and—unlike here—the isotherms in earlier studies were often shifted artificially to match the experimental ones (Baoukina et al., 2007b; Huynh et al., 2014). This signals that, in previous simulation studies of pulmonary surfactant monolayers, the descriptions of the physics at the interface was inadequate due to inadequate water models (Baoukina et al., 2007b; Lamberg and Ollila, 2015). The simulation approach used in the present work largely overcomes these issues through a combination of recently developed simulation models and simulation parameters (Klauda et al., 2010; Izadi et al., 2014; Javanainen et al., 2017a). It is worth noting that the majority of atomistic water models also feature significant errors in air–water surface tension (Chen and Smith, 2007; Vega and De Miguel, 2007), and that OPC provides not only an accurate description of the surface tension, but also of many key properties of water (Izadi et al., 2014). As such performance is not achievable using three-site models (Izadi and Onufriev, 2016), and since the 4-point OPC water does not compromise the behavior of CHARMM36 lipids (Javanainen et al., 2017a), this combination should perhaps be adopted more widely.

The only apparent deviations between the current simulation model and experiment appear at small APLs, where a monolayer collapse plateau is observed in the experimental isotherms at above 45 mN/m and 50 mN/m surface pressures (with areas below 52 and 54 Ų) at 298 and 310 K, respectively. This equilibrium collapse pressure is in the ballpark measured for numerous single-component monolayers (Mansour and Zografi, 2007). In the simulations at 298 and 310 K, the systems remain stable above the experimental equilibrium collapse pressures of phospholipid monolayers (Lee et al., 2001; Mansour and Zografi, 2007), which likely results from the fact that the simulated monolayer is kinetically trapped in a metastable state and its collapse is limited by the system size and its periodic nature. This limitation has two implications. On the one hand, it hinders studies of the formation of bilayer folds upon high lateral compression, i.e., during collapse. For such studies, the qualitative picture provided by the MARTINI model (Marrink et al., 2007) is likely sufficient (Baoukina et al., 2007a, 2008, 2014). On the other hand, the atomistic simulations might model well the non-equilibrium situation in the lungs, where monolayer collapse is prevented by rapid compression, and surface pressures can thus reach very high values (Crane and Hall, 2001). Experimentally, achieving such metastable states requires a specific Langmuir trough or the use of captive bubble surfactometer (Schürch et al., 1998; Crane et al., 1999).

One obvious question that arises is whether our AFM observations of sub-micrometer-sized domains are compatible with the nanoscopic and fairly transient domains detected in our simulations. This is a general problem in atomistic simulations, including those dealing with lipid bilayers that often compare heterogeneity between a simulation system with a few hundred lipids to experimental data obtained on phase-separating systems, where the phases would coalesce in any size scale if given enough time (Sodt et al., 2014). For us the situation is slightly more favorable, as the length scales in our simulations and experiments are only 1–3 orders of magnitude apart. Line tension, i.e., the penalty of creating a domain boundary, dictates domain shapes (Bernardino de la Serna et al., 2004). The non-circular domain shapes observed by AFM indicate a fairly small line tension, which is also supported by the fact that small and independent domains do not coalesce—as would eventually happen in the case of phase separation. This is in line with the short lifetimes of domains observed in our simulations (Figure 5C). Based on the simulations, we hypothesize that the L_C-like domains are highly dynamic in monolayers at the air–water interface, but much less dynamic when the samples are transferred onto a mica substrate for AFM measurements. Nevertheless, the heterogeneity observed by AFM is not an artifact of the immobilization, as heterogeneity was also detected in the dilatational modulus of monolayers measured in the Langmuir trough (Figure 3C). Finally, it is worth again noting that the domains visible in our simulations are limited by the simulation box size, despite representing the state-of-the-art in this regard. Our control simulations with more lipids (see section 2) result in similar behavior, thus suggesting that the observed heterogeneity does not arise due to finite-size effects. The surface pressures were also unchanged by the change in simulation box size (Supplementary Figure 1). Still, the domain morphologies are likely dependent on the system size, and unfortunately proper lateral mixing of lipids is not feasible in the larger monolayers with atomistic simulations, nor in the smaller monolayer when most of the monolayer is part of the L_C-like regions. For the two key functions—the viscoelastic properties of the surfactant and phase-dependent protein functions (Bernardino de la Serna et al., 2004, 2009, 2013a,b; Pérez-Gil, 2008; Casals and Cañadas, 2012)—the domain morphology is likely not very important. The domain lifetimes extracted from simulations are likely dependent on the domain sizes, and thus also system sizes, since escaping a core of a larger domain takes more time. Unfortunately, extrapolating these lifetimes to microscopic scales is not straightforward.

Our study highlights that, despite the absence of a coexistence plateau in the surface pressure-area isotherm, heterogeneity at the nanometer scale—undetectable by fluorescence or Brewster angle microscopies—should be detected by experimental approaches other than AFM. The surface dilatational rheology measurements using an oscillating barrier readily observed a major change in the dilatational modulus at an APL where the heterogeneities begin to appear (Figure 3C). Lipid diffusion analyzed from the simulations also revealed a change in the trend at this crossover APL. As monolayer phase transitions are readily detected in experiments that probe lipid diffusion (Peters and Beck, 1983), such experiments could also detect smaller scale heterogeneity. Finally, lipid tilt measured from our simulations (Figure 3B) revealed a persistent tilt angle of 25° up to the crossover APL, and this behavior should be captured by either X-ray diffraction (Lee et al., 2002) or vibrational sum frequency generation spectroscopy methods (Ma and Allen, 2007).

5. Conclusions

Using a combination of Langmuir trough experiments, AFM imaging, and atomistic molecular dynamics simulations, we demonstrated that a synthetic quaternary lipid mixture is able to qualitatively reproduce the key features of the phase behavior of the native pulmonary surfactant extracts. Under a large range of compression levels, thicker L_C-like domains appear in the otherwise thinner L_E-phase monolayer. These domains are dynamic and only slightly enriched in DPPC with two saturated chains. We demonstrated that despite there not being a visible phase transition to the L_C phase, some monolayer properties change significantly at well-defined values of area per lipid, and this crossover value is consistent between numerous quantities. Moreover, since these properties should be readily measurable using experimental methods, our study also guides experimental work on detecting heterogeneities in biofilms.

Synthetic pulmonary surfactants mimicking the properties of the full functional surfactant are continuously investigated (Walther et al., 2019). A critical aspect is to find a proper lipid composition mixture. Our results highlight how already a few key lipid components of the pulmonary surfactant display small domains, resembling the behavior of surfactant extracts (Bernardino de la Serna et al., 2013a). The lipid mixture is able to pack in a dynamic manner, thus enabling efficient surface tension reduction while maintaining sufficient fluidity. This behavior might also be crucial for the function of surfactant proteins, as has been investigated in other molecular simulations with multilipid components (Robichaud et al., 2019), which we will focus on in our future work.

Our approach, combining the CHARMM36 force field with the four-point OPC water model, enables atomistic studies of lipid structures at the air–water interfaces in the complex pulmonary surfactant, allowing for studies of the physiologically important processes in the lung at a detail difficult to achieve experimentally. By integrating experimental data with molecular simulations, we provide, for the first time, a quantitatively accurate, unprecedented picture of the structural and dynamic properties of a realistic model of lung surfactant, under physiologically relevant conditions.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Author Contributions

MJ, LM, and IV designed the simulations. JL performed all simulations under the supervision of MJ and IV. JL and MJ analyzed the simulations. JB designed the isotherm and AFM experiments and analyzed the data. BS performed these experiments under the supervision of JB. RP designed, performed, and analyzed the surface dilatational rheology experiments. JL and MJ wrote the first draft of the manuscript. All authors contributed to manuscript revision, read, and approved the submitted version.

Funding

JL and IV thank the Sigrid Juselius Foundation and the Academy of Finland (Centre of Excellence program) for financial support. IV also thanks the HiLIFE (Helsinki Institute of Life Science) Fellow program. RP thanks the Mary and Georg C. Ehrnrooth Foundation, the Eye and Tissue Bank Foundation, the Evald and Hilda Nissi Foundation and the Biomedicum Helsinki Foundation for financial support. LM acknowledges funding by the Institut national de la santé et de la recherche médicale (INSERM). MJ thanks the Emil Aaltonen Foundation for funding. BS and JB acknowledge funding from the Bill and Melinda Gates Foundation, Grant No. INV-016631.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank CSC–IT Center for Science (Espoo, Finland) for computing resources.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2020.581016/full#supplementary-material

References