Domain size polydispersity effects on the structural and dynamical properties in lipid monolayers with phase coexistence

Lipid monolayers with phase coexistence are a frequently used model for lipid membranes. In these systems, domains of the liquid-condensed phase always present size polydispersity. However, very few theoretical works consider size distribution effects on the monolayer properties. Because of the difference in surface densities, domains have excess dipolar density with respect to the surrounding liquid expanded phase, originating a dipolar inter-domain interaction. This interaction depends on the domain area, and hence the presence of a domain size distribution is associated with interaction polydispersity. Inter-domain interactions are fundamental to understanding the structure and dynamics of the monolayer. For this reason, it is expected that polydispersity significantly alters monolayer properties. By means of Brownian dynamics simulations, we study the radial distribution function (RDF), the average mean square displacement and the average time-dependent self-diffusion coefficient, D(t), of lipid monolayers with normal distributed size domains. It was found that polydispersity strongly affects the value of the interaction strength obtained, which is greatly underestimated if polydispersity is not considered. However, within a certain range of parameters, the RDF obtained from a polydisperse model can be well approximated by that of a monodisperse model, suitably fitting the interaction strength, even for 40% polydispersities. For small interaction strengths or small polydispersities, the polydisperse systems obtained from fitting the experimental RDF have an average mean square displacement and D(t) in good agreement with that of the monodisperse system.


Introduction
Most biologically relevant monolayers present phase coexistence characterized by domains formed by lipids in an ordered phase state, dispersed in a continuous, disordered phase [1][2][3][4] . The domains interact with each other [5][6][7] , and these interactions affect their own movement 8,9 as well as that of other species present in the monolayer 10,11 . Inter-domain interaction may be related to electrostatic forces (dipolar or Coulombic repulsions), forces related to the spontaneous curvature of the coexisting phases and hydrodynamic forces that appear when domains are in motion.
Dipolar inter-domain interaction, which is always present, arises from the excess dipolar density of the ordered phase with respect to the continuous phase. Hence, the dipolar strength is proportional to the domain area.
Lipid monolayers typically exhibit domain size polydispersity [12][13][14][15] . In particular, Langmuir monolayers at the air-water interface show a wide domain size distribution 9,10,14,16 . Due to the origin of the dipolar interaction, size polydispersity leads to interaction polydispersity, which usually turns out to be the most important, as a consequence of the quadratic dependence of dipolar strength on domain size.
Structural and dynamical properties of monolayers are mainly determined by inter-domain interactions and hence, for systems where the dipolar interaction is dominant, are expected to be strongly affected by the presence of domains of different sizes.
To the best of our knowledge, the effects of size polydispersity on the structural and dynamical properties of lipid monolayers have not been studied. However, for the determination of physical parameters of the constituting lipids, size distribution has been taken into account. Mulder 12  cular domains in Langmuir monolayers for determining physical parameters of surfactants. He approximates the exact size distribution by a Gaussian distribution and uses a simplified theoretical analysis, in which the inter-domain interactions are treated approximately by considering equally sized domains arranged in a regular hexagonal array. Lee et al. 14 obtain the excess dipolar density by fitting the size distribution with an equilibrium thermodynamic expression. Their scheme assumes no interactions between domains, and hence it is valid for sufficiently diluted systems.
In a previous work 16 , we proposed a novel way of estimating dipolar repulsion, using a passive method that involves the analysis of images of the monolayer with phase coexistence. The method is based on comparing the pair correlation function obtained from experiments with that obtained from simulations of systems of monodisperse domains interacting by a dipolar density pair potential. We also studied the point dipole approximation for the dipolar density pair potential, and determined an effective point dipole interaction strength that reproduces the structural properties of a system with dipolar density pair potential.
In this work, we use Brownian dynamics simulations to investigate the effects of polydispersity on the structure and dynamics of lipid monolayer models, with parameters chosen within a range of experimental interest. For this purpose, we study the pair cor-relation function, the average mean square displacement and the average time-dependent self-diffusion coefficient, as a function of polydispersity and dipolar interaction strength. In addition, we study how polydispersity would affect determination of dipolar repulsion strength from experimental data, when the method proposed in Ref. 16 is used based on a monodisperse model.

Polydisperse domain interactions
We consider a monolayer in its two-phase liquid-condensed (LC) and liquid-expanded (LE) coexistence region, where the LC phase forms domains in the LE phase, which occupies the larger area of the monolayer. Because of the difference in surface densities, the LC domains possess an excess dipole density, σ , with respect to the surrounding LE phase 42 . This originates dipolar repulsive interactions between the domains. We model the mixed monolayer as a 2D-dispersion of domains which interact through a dipolar pair potential. Considering only dipole components perpendicularly oriented to the interface, and using the approximation of a point dipole in the center of each domain with an excluded area, the resulting pair potential between domain i and domain j can be described by: where u hc (r) is a hard core repulsive potential, and where µ i is the dipole moment of the domain i representing the dipole density σ over the domain area A i , µ i = σ A i , r center-tocenter domain distance, ε 0 is the vacuum permittivity and ε * is an effective permittivity 43 that considers the relative permittivities of the membrane, ε m , water, ε w , and air, ε a : We consider polydisperse circular domains of radius R i , accordingly Eq. 2 becomes: This equation makes explicit how domain size distribution leads to interaction polydispersity. For convenience, we define a dimensionless interaction strength: where k B is the Boltzmann's constant, T the absolute temperature and R m the mean value of the radii distribution. Then, the dimensionless dipolar pair potential takes the form:

Domain size distribution
Most experiments on lipid monolayers present domain size polydispersity 4,8,13,14,16 . The functional form of the domain size distribution is clearly crucial to the structural and dynamical properties of the monolayer. Here, we approximate the size polydispersity by a truncated normal distribution function, where R m is the mean or expectation of the distribution, Σ is the standard deviation, and a is a normalization constant which arises from simetrically truncating the normal distribution function to exclude negative domain radii. The system polydispersity, ω, is characterized by the ratio of the width of the distribution to its mean; ω = Σ/R m . In our simulations we use a discrete counterpart of P(R) (see Fig. 1), where we chose a number of domain species N pol , and a bin width s, such that the radii distribution of the system is described by the set {R α , N α /N}, α = 1, N pol , where N α is the number of domains of type α and N is the total number of domains. For the system with the largest polydispersity studied, ω = 0.4, 1% of the domain radii fall outside the histograms.

Radial distribution function
A key quantity for characterizing the structure of the monolayer is the radial distribution function (RDF), g(r). Considering a distribution of domains in the monolayer plane, g(r) is related to the probability of finding a domain at a distance r from another domain chosen as a reference point: where g α,β (r) are the partial radial distribution functions, defined as: with A = L 2 the total monolayer area and the angular brackets indicating an equilibrium ensemble average. Note that ρ β g α,β (r) is the probability density of finding a β particle at a distance r from an α particle, where ρ β = N β /A is the number density of domains with radius β .

Diffusion
In order to evaluate the effects of polydispersity on the system dynamics, we studied the mean square displacement (MSD) and selfdiffusion of domains. The mean square displacement of a domain of type α with its center at position r 1,α (t) at time t is given by: where the angular brackets indicate an equilibrium ensemble average. However, many experimental results do not differentiate between domain radii. Hence, we evaluate the average MSD as a representative quantity: where x α = N α /N is the molar fraction of domains of type α. The time-dependent self-diffusion coefficient, D α (t), is defined as the time derivative of W α (t). The short-time limit of D α (t) corresponds to the free domain diffusion coefficient, which in this work is approximated by the diffusion coefficient of a disk in a two-dimensional simple fluid: where η is the viscosity of the fluid. It is important to note that, in our simulations, only the ratios between the different D 0 α are relevant. The long-time self-diffusion coefficient is defined as: Analogously, the average time-dependent self-diffusion coefficient, D(t), is obtained from the time derivative of W (t), and the corresponding limits are:

Simulations
We model the mixed monolayer as a two-dimensional Brownian suspension of interacting hard disks with polydisperse radii immersed in an effective fluid, each disk representing an idealized lipid domain. The inter-domain interactions are described by the point dipole pair interaction Eq. (6) plus a hard disk repulsive part. Hydrodynamics interactions are disregarded. To study the static and dynamical properties of the mixed monolayer model, we performed Brownian dynamics (BD) simulations. In this scheme, the finite difference equation describing the in-plane displacement of N Brownian disks immersed in a fluid during the time step ∆t is given by 44 where r α,i (t) is the position of domain i of type α at time t, F P β , j is the force on disk i of type α due to the disk j of type β and X α,i a random displacement vector of domain i of type α originating from solvent particle collisions. X α,i is sampled from a Gaussian distribution with zero mean and covariance matrix: where I is the identity matrix, and δ i, j the Kronecker delta. The simulated systems consisted of N disks with radii distribution {R α , N α /N} under periodic boundary conditions, using the minimum image convention.
The size of the simulation box, L, was determined from the condensed area fraction, φ , defined as: The system is completely characterized by the following parameters: the interaction strength f (or equivalently, the dipolar density σ ), the size distribution {R α , N α /N} and the total area fraction φ .
Throughout this work, the area fraction is fixed at the value φ = 0.20, which was chosen as a typical value of experimental monolayer micrographs.
In Figure 1, we show the distributions used in our studies and, in Table S1 of the electronic supplementary information (ESI), the size distribution {R α , N α /N} of each system studied can be found. We have carried out simulation studies with four distinct polydispersities: ω = 0.1, ω = 0.2, ω = 0.3 and ω = 0.4. These systems are compared with simulations of perfectly monodisperse domains.
The bin width was chosen as s/R m = 0.2. This choice leads to a small number of domain types for each polydispersity, and at the same time reproduces the distribution shape qualitatively well. For a typical experimental average radius of R m = 5 px, the bin width s = 1 px is in the range of typical optical microscopy experimental error. We verified that, for the systems studied, the Gaussian distribution is well-described by the selected discretization {R α , N α /N}. The RDF and the average MSD do not change substantially if half of the bin width value is used.
Finally, the interaction strength f was varied such that the system remains in its fluid phase and has experimental interest.
The  Table S1 of the ESI. For the systems studied, we verified that there is no system size dependency.

Effects of polydispersity on g(r) and MSD(t)
With the aim of analyzing the effect of polydispersity, we consider systems with four different interaction strengths, which correspond to the liquid regime for the selected area fraction (φ = 0.2), and for each interaction strength we vary the polydispersity from ω = 0 (monodisperse) to 0.4. Here, we present the results for f = 1.2 and 4.8, and, in the ESI, the intermediate values of f = 2.4 and 3.6 are shown (Fig. S1). Figure 2 shows the pair correlation function for f = 1.2 (a) and 4.8 (b), for ω = 0, 0.1, 0.2, 0.3 and 0.4. In both cases, we observe a decrease and a broadening of the first peak as the polydispersity grows. Furthermore, there is a shift of the peak position to greater distances and the g(r) start to show correlations for shorter distances. In general, polydispersity softens the peaks and minima of the g(r). It is remarkable that, even for %40 of polydispersity, there is still a well-developed first minimum and second maximum, for highly interacting systems. This can be attributed to the fact that the interaction strength grows with the fourth power of the domain size.
To analyze how polydispersity affects the dynamical properties of the monolayer, we calculated the average MSD and the average time-dependent diffusion coefficient. Results for the selected systems are presented in Fig. 3. For the weak interaction f = 1. ( Fig. 3a), we observe that the average short-time diffusion coefficients of polydisperse systems are larger than that of monodisperse system. This is a direct consequence of the inverse radius dependence of the short-time coefficient diffusion, D 0 α , and of the symmetry of the radii distribution. Here, for the more polydisperse case, ω = 0.4, the difference is about 20%. For the average long-time diffusion coefficients we find a similar ordering, but with slightly larger differences. In particular, for ω = 0.4 the difference reaches 40%. However, the relative decrease of the average long-time diffusion constant with respect to the short-time limit, D L /D 0 , differs by less than 15%. For all polydispersities, D(t)/D 0 diminish roughly 50% at the long-time limit. In the inset of Fig. 3a, the respective average MSDs are shown. Note that, within the logarithmic scale, the polydisperse systems appear very similar to the monodisperse system, except for a slight translation to larger values. Figure 3b shows the average diffusion quantities for f = 4.8. For these systems, the same qualitative behavior is observed as in previous ones. In particular, the short-time limits are identical since they have identical radii distribution. In addition, the average long-time diffusion coefficients are smaller, because of the stronger interactions.
Comparing the most polydisperse system, ω = 0.4, with the monodisperse one, it is found that the average long-time diffusion coefficient is of the order of 60% larger, while the ratio D L /D 0 is roughly 30% larger. In general, an enhancement of average diffusion is observed as polydispersity increases. This effect is stronger in the long-time regime and is also more pronounced for systems with stronger interactions.
To describe the structure of the polydisperse system in more de- tail, we show, for the system with ω = 0.3 and f = 4.8, the partial radial distribution functions g α,α (r) in Fig. 4a and g 1,α (r) in Fig. 4b. These show that the partial RDFs start to differentiate from zero at distances larger than the respective contact values, i.e., domains do not come into contact. Note that the distance between the first and the second neighbor shell and the first minimum depth are very similar for different domain sizes. This is clearly seen in the inset of Fig. 4a, where we have plotted g α,α (r) horizontally shifted by their respective peak positions, r max . Considering the spatial correlation between the smaller domains and the other domain types, as shown in Fig. 4b, it is observed that the probability of finding the smallest domains (type α = 1) around domains of other types increases with domain sizes. However, the first minimum depth of g 1,α (r) decreases with domain type α. The other g α,β (r), not shown here, also indicates that small domains are more probably found close to larger ones. Focusing now on the time dependent diffusivities, in Fig. 5 we show the time dependent diffusion coefficient for each domain type, D α (t), and the corresponding average, D(t). Here, to compare the relative slowdown of the dynamics between the different types, the diffusivities are normalized by their corresponding short time limit, D 0 α . It is observed that, except for the smallest domain type, the self-diffusion of the different domain sizes slows down similarly, reaching a long-time limit of roughly D L α /D 0 α ≃ 0.29. There is a remarkably large difference in the behavior of the smallest domains, which show a much smaller slowdown than the others, i.e., they are able to leave the neighbor cages more easily. This could be attributed to the fact that the dipolar strength of the do- mains grows quadratically with the radius.

Monodisperse vs. polydisperse models
In lipid monolayers with phase coexistence, the dipolar repulsion σ ( f ) is usually not a known parameter. One method to determine it is to fit the experimentally measured RDF with one obtained from simulations, using σ as the only adjustable parameter 16 . This assumes monodisperse distribution of domain radii, and hence monodisperse inter-domain interactions. Therefore, it is important to assess the effects of polydispersity on the determination of σ using this method. For this purpose, we used the g(r) from the monodisperse systems studied in this work as the reference RDF, as if they were previously fitted to experimental data sets. Then, we fitted the reference g(r) with the polydisperse systems shown in Fig. 1 (as was already stated φ = 0.2 for all systems) and we analyzed how the interaction strength varies with ω. Note that, as in our previous work 16 , only the first peak height of the reference g(r) is considered in the fitting procedure. Figure 6 shows the results for systems with g(r max ) = 1.35 (a) and g(r max ) = 2.01 (b) (monodisperse models with f m = 1.2 and 4.8, respectively). In the ESI, systems with g(r max ) = 1.63 and g(r max ) = 1.84 ( f m = 2.4 and 3.6, respectively) are shown (Fig. S3).
The RDF from the monodisperse system with g(r max ) = 1.35 (Fig. 6a) is qualitatively well captured by the polydisperse systems, except for short distances, where for more polydisperse systems the g(r) start to grow at shorter distances, as expected. In the inset, we have plotted the g(r) horizontally shifted so that the first peak position coincides. This clearly shows that the distance between the different neighbor shells and the depth of minima are very similar for all polydispersities. This agreement is striking, since we are comparing systems with polydispersities as large as 40%.
For a more structured system with g(r max ) = 2.01, shown in Figure 6b, we observe similar behavior. However, for polydispersities larger than 20%, the depth of the first minimum starts to differentiate, at the same time as the third peak position begins to dephase. This can be seen in the inset. Note that, for 40% polydispersity, the fit is already not possible, i.e., there is no interaction strength for which the resulting g(r) reaches the maximum value of 2.01.
In general, to reach a certain value of g(r max ) a larger interaction strength is needed as the polydispersities increase. However, not all the values of g(r max ) can be obtained for systems with large polydispersity. Figure 7 shows g(r max ) as a function of the interaction strength f for the different polydispersities studied here. It can be seen that the values g(r max ) tend to saturate for high interaction strength.
For experimental systems with φ = 0.20 and approximately Gaussian size distributions, Figure 7 can be used as a working curve to estimate f (or σ ) directly from the experimental g(r) without implementing any simulation, generalizing the method presented in Ref. 16 to polydisperse systems.
At this point, we introduce a new parameter Γ, which combines the number density, ρ = N/A, and the interaction strength, f , in one independent dimensionless parameter. Namely, where r m = ρ 1/2 is the mean geometrical distance between domains. Systems for which the hard disk interactions can be disregarded (i.e., strongly interacting or low density systems) are completely determined by Γ, the scaled domain radii of the different species, λ α = R α /R m , and their corresponding molar fractions, Given that, for the systems studied here, the hard disk interaction turned out to be irrelevant, Fig. 7 may be recast to show the maximum of the RDF as a function of the Γ, as shown in Fig. 8. In this way, if the size polydispersity can be described by a Gaussian distribution with {ω, R m }, and g(r max ) and ρ are determined from the experiments, the Γ value could be estimated directly from this figure, and subsequently the dipolar density can be obtained.
To further study how polydispersity affects the structure, Figure 9 shows the interaction strength needed to obtain a certain g(r max ) as a function of ω. The curves for g(r max ) = 1.35 and g(r max ) = 2.01 correspond to the RDFs shown in Fig. 6 a and b, respectively. Here, it can be clearly seen that the interaction strength for a given g(r max ) increases notably faster for more structured systems and, in particular, for g(r max ) = 1.35 and ω = 0.3, f / f m = 2 and for g(r max ) = 2.01 and ω = 0.3 f / f m = 5.21. The dotted line for the case g(r max ) = 2.01 indicates that no system with ω = 0.4 is able to reach this peak height.
Finally, we consider the dynamics of the systems studied shown in Figure 6. The corresponding average diffusion coefficient and MSD are shown in Figure 10.
It is remarkable that, for the less structured system (Fig. 10a), polydispersity does not much affect the intermediate and longtime average dynamical quantities. However, a completely different scenario occurs for the more structured system (Fig. 10b). In this case, up to ω = 0.1, the average diffusion coefficient behaves similarly to that of the monodisperse system. On the other hand, for larger polydispersities, D(t) strongly deviates from the monodisperse system. The system reaches the sub-diffusive regime faster, as polydispersity increases. Besides, the long-time average diffusion coefficient is smaller for larger polydispersities. This is probably a consequence of the dependence of the interaction strength on the domain sizes, which is evident for strongly interacting systems.

Conclusions
In this work, we studied the influence of domain-size polydispersity on the structure and dynamics of model lipid monolayers with phase coexistence at the air-water interface. The sizepolydispersity was modeled as a discretized Gaussian distribution. Studying the monolayer structure, we found a decrease and a broadening of the first peak of the RDF as the polydispersity grows. Furthermore, a shift of the peak position to greater distances and the occurrence of correlations for shorter distances were observed. Notably, highly interacting systems ( f = 4.8) presented a welldeveloped first minimum and second maximum, even for 40% of polydispersity. For all the systems studied, the partial RDF starts to differentiate from zero at distances larger than the respective contact values. This indicates that domains do not come into contact and hence that the hard-disk interaction is not relevant. Regarding the spatial correlation between the smaller domains and the other domain types, it was observed that the probability of finding the smallest domains around domains of other types increases with domain sizes.
Analyzing the domain dynamics, we found an enhancement of the average diffusion as polydispersity increases. This is more pronounced in the long-time regime. For systems with stronger interactions, the overall enhancement is more noticeable. It was also found that the self-diffusion of the smallest domains shows a much smaller slowdown than the other domain sizes, i.e., they are able to leave the neighbor cages more easily. Interaction strength of the polydisperse models that lead to a similar structure to that of the monodisperse model, for different monodisperse systems. The lines are a guide for the eye.
We also studied the effects of polydispersity in the determination of f (or σ ) by the method proposed in Ref. 16 , where a monodisperse model is used to fit the experimental RDF. For the systems considered, it was found that polydispersity strongly affects the value of f obtained, which is greatly underestimated if polydispersity is not considered. Only for the experimental system with small polydispersities and/or weak interactions is achieved a good approximation.
It is remarkable that, even for large polydispersities, the fitted RDFs result in good agreement with the reference one; they have the same second and third peak heights and neighbor shell distances. They mainly differ in the depth of the first minimum and for distances where the g(r) starts to grow.
With regard to the dynamics, on the other hand, only for weak interactions or small polydispersities does the average timedependent diffusion coefficient agree with the reference system. For stronger interactions or larger polydispersities, a noticeable slow-down is observed in the average dynamics.
Finally, the method proposed in Ref. 16 may be straightforwardly generalized to include polydispersity by directly fitting the experimental g(r) with a model that accounts for the measured domain size distribution in the simulations. Alternatively, using a Gaussian size-distribution model and for a selected range of { f , ω}, a set of figures, like Fig. 7, can be generated for different values of φ , for later use as working curves to estimate f (or σ ) directly from the experimental data without implementing any simulation. In particular, for a system in which the hard disk interaction can be disregarded, only one working curve is needed (Fig. 8).