Molecular simulations of micellar aggregation of polysorbate 20 ester fractions and their interaction with N-phenyl-1-naphthylamine dye

    Molecular simulations of micellar aggregation of polysorbate 20 ester fractions and their interaction with N-phenyl-1-naphthylamine dye Mauro Lapelosa, Thomas W. Patapoff, Isidro E. Zarraga PII: DOI: Reference:

S0301-4622(16)30049-7 doi: 10.1016/j.bpc.2016.03.003 BIOCHE 5888

To appear in:

Biophysical Chemistry

Received date: Revised date: Accepted date:

26 February 2016 22 March 2016 22 March 2016

Please cite this article as: Mauro Lapelosa, Thomas W. Patapoff, Isidro E. Zarraga, Molecular simulations of micellar aggregation of polysorbate 20 ester fractions and their interaction with N-phenyl-1-naphthylamine dye, Biophysical Chemistry (2016), doi: 10.1016/j.bpc.2016.03.003

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Molecular simulations of micellar aggregation of polysorbate

SC

RI

naphthylamine dye

PT

20 ester fractions and their interaction with N-phenyl-1-

NU

Running title: polysorbate 20 micellar aggregation

MA

Late Stage Pharmaceutical Developmenta, Early Stage Pharmaceutical Developmentb, Genentech Inc. (member of Roche), South San Francisco, California 94080, USA

AC CE P

TE

D

Mauro Lapelosaa*, Thomas W. Patapoffb, and Isidro E. Zarragaa

* Corresponding author: Late Stage Pharmaceutical development, Genentech Inc. (member of Roche), 500 Forbes Blvd, South San Francisco, California 94080, USA E-mail addresses: [email protected] (M. Lapelosa)

1

PT

ACCEPTED MANUSCRIPT

Abstract

RI

Micellar aggregation behavior of polysorbate 20 (PS20) has generated significant interest

SC

because of the wide use of PS20 as a surfactant to minimize protein surface adsorption

NU

and mitigate protein aggregation. Thus, there is a need for better molecular understanding of what drives the biophysical behavior of PS20 in solution. We observe that a complex

MA

amphipathic PS20 molecule, which contains both hydrophobic tail and relatively large hydrophilic head, self-associates strongly within the course of a molecular dynamics

D

simulation performed with a fully atomistic representation of the molecule and an explicit

TE

water solvent model. The in silico behavior is consistent with micellar models of PS20 in

AC CE P

solution. The dynamics of this self-association is rather complex involving both internal reorganization of the molecule and diffusion to form stable micelle-like aggregates. The micellar aggregates of PS20 are long-lived and are formed by the balance between the large hydrophobic interactions associated with the aliphatic tail of PS20, and the steric repulsion of the hydrophilic sorbitan head structure. In the present work, molecular models of PS20 that represent naturally occurring PS20 fractions were produced and characterized in silico. The study investigated the monoester and diester fractions: PS20M, and PS20D. These fractions present differences in the strength of their hydrophobic effect, which influences the aggregation behavior. Adaptive biasing force (ABF) simulations were carried out with the PS20M and PS20D molecular constructs to calculate the free energy of their pairwise interaction. The free energy barrier for the dissociation is higher for PS20D compared with PS20M. The results show that hydrogen

2

ACCEPTED MANUSCRIPT bonds can form when head groups are in close proximity, such as in the PS20 aggregate assembly, and the free energy of interaction can be used to predict the morphology of the

PT

micellar aggregate for the different PS20 fractions. We were also able to simulate PS20

RI

in the presence of N-phenyl-1-naphthylamine (NPN) to study the solution behavior of the

NU

hydrophobic core of the PS20 micellar aggregate.

SC

hydrophobic molecule and of the mechanism in which it is sequestered in the

MA

Keywords: polysorbates, molecular dynamics, polysorbates 20 aggregation, hydrophobic

AC CE P

Introduction

TE

D

effect

Polysorbates are widely used in protein formulations to enhance protein stability. They undergo one of the most common phenomenon occurring in surfactants: micellar aggregation [1,2]. Polysorbate 20 (PS20) is comprised of a polyoxyethylene sorbitan head and an aliphatic tail. Ethylene oxide chains are attached at three different hydroxyl positions in the sorbitan moiety. The number of repeat ethylene oxide unit can vary in the specific position, but the total number sums up to 20 (x +y +z +w) [1]. The PS20 is a mixture of various esters, but the largest fraction (~40-60 %) is comprised of lauric acid esters. For simplicity, our simulations focus here on PS20 fractions associated with the most prominent fixed-length aliphatic chain, i.e. the monolaurate and di-laurate esters. For a monolaurate, the fatty acid moiety is linked through an ester bond to the ethylene

3

ACCEPTED MANUSCRIPT oxide oxygen at the z position. The particularly strong amphipathic properties of PS20, having a very hydrophilic head and very hydrophobic tail, make it for a great surfactant

PT

to use for a wide range of protein formulation applications, e.g. to prevent protein

RI

aggregation and adsorption on interfaces (e.g. air-liquid, oil-water). The ability to mitigate protein aggregation is important, for example, in ensuring the quality of a

SC

biopharmaceutical product despite exposure to interfacial stresses during its production,

NU

shipment, and storage.

The existence of micellar aggregates is inferred from surface tension measurements, but

MA

such measurements lack structural understanding at the molecular level. The amphipathic nature of PS20 molecules, together with the relatively large hydrophilic head group for

TE

D

example, can cause the hydrophobic tail to arrange so that exposure to the solvent is minimal (e.g. with hydrophilic moiety surrounding the tail, or with the help of other PS20

AC CE P

molecules in a micelle-like aggregate). The actual molecular details of such individual PS20 molecule conformation are not known. It is therefore important to use molecular dynamics for characterizing the conformational ensemble of PS20 as well as for evaluating the hydrophobic effect in the micelle aggregation process. Recently, however, the large availability of computational power has enabled atomistic molecular dynamics (MD) simulations to be carried out for systems containing up to millions of atoms allowing to study aggregation in fully atomic details over time scales of hundreds of nanoseconds [3]. MD simulations provide an atomic resolution of structural details down to the Angstrom level, and reveal dynamics and possible kinetics of molecular processes relevant to aggregation of surfactants in water [4]. This provides the motivation to use atomistic MD simulations to study aggregation properties of surfactants, as well as

4

ACCEPTED MANUSCRIPT surfactant behavior at air-liquid interfaces, and for understanding the structural details within a micellar aggregate.

PT

All-atom MD simulations based on effective molecular mechanics potentials have been

RI

used to study the micelles of octyl glucoside [5], C12E4 [6], dodecylphosphocholine (DPC) [4,7,8], polysorbate 80 [9], sodium dodecyl sulfate (SDS) [10,11,12] and sodium

SC

octanoate [13,14]. More sophisticated free energy methods can measure the free energy

NU

of interaction [15,16,17,18,19,20,21,22,23,24,25,26], in particular adaptive biasing force (ABF) simulations [27] have been proved successful in binding problems such as

MA

association of transmembrane α-helices [28]. The method approximates the local mean force at ξ(q) = z by averaging over instantaneous forces observed by the system

TE

D

trajectory. This estimated mean force, A′(t), is then subtracted from the dynamics to reduce and eventually eliminate any force along a collective variable (ξ) [27,29]. This

AC CE P

method turns out to be well suited to measure the free energy of interaction for the PS20 molecules with sorbitan monoester and sorbitan diester tails. The hydrophobic effect is the phenomenon of solvent induced interaction between two or more apolar molecules [30,31,32,33,34,35]. These interactions are the main driving forces responsible of protein folding, biopolymer interaction, stability of micelles and membranes in aqueous solutions. The formation of clusters for hydrophobic particles in water is explained in terms of the dependence of solvation structure on solute size. If we place a number of small particles in water, separated by a large distance from each other, then the solvation free energy will simply be the sum of the free energy contribution for each solvent particle. However, when the solute molecules aggregate together to form a cluster or micelle [30], and the free energy grows linearly with exposed surface area, this

5

ACCEPTED MANUSCRIPT suggests the clustering is driven by hydrophobic interactions. In contrast, electrostatic interactions are more dependent on ionic strength, net charge, and electric multipole

PT

effects. The solvation free energy for hydrophobically interacting molecules reaches a

RI

lower value than the sum of solvation free energies for each molecule as self-association progresses. Well documented studies examining the standard free energy of transfer for

SC

many hydrocarbons compounds from gas phase to water and from gas phase to a

NU

hydrophobic solvent provided evidence that the number of hydrogen-carbon bonds represents the basis for the hydrophobic effect [34].

MA

The aggregation of polysorbates to form a micelle, as mentioned earlier, is related to the hydrophobic effect. The general model inferred from surface tension experiments shows

TE

D

that PS20 molecules form a monolayer at the air/water while minimizing exposure of the hydrophobic tail to water. These can occur at low concentrations of surfactant, however

AC CE P

as the critical micelle concentration (CMC) [36] is approached the amount of surfactant the air-liquid surface can accommodate saturates and the surface tension plateau provides evidence of this. In the bulk, such a critical concentration implies the surfactant molecules are at close enough proximity to spontaneously phase separate into a micellar aggregate. Thermodynamically, the related chemical potential allows the surfactant molecules to organize themselves in solution with the hydrophilic groups oriented towards the bulk water and the hydrocarbon chains pointing towards the interior of the hydrophobic core with minimal or no water [37]. Molecular understanding of this phenomena have been attempted in previous MD studies on some ionic surfactants, but no extensive study has been conducted for nonionic polysorbate, especially its associated ester fractions. In previous studies, MD simulations

6

ACCEPTED MANUSCRIPT of sodium dodecyl sulfate (SDS) were able to show the formation of a complete micelle in water [12] and around a protein.

The stability of this protein in the micelle

PT

environment was observed from an initial random placement of SDS molecules in an

RI

aqueous environment [38].

This study reports the use of molecular dynamics simulations with fully atomistic details

SC

with explicit solvent using the recent version of the CHARMM force field [39] to study

NU

the aggregation of PS20 into micelle units from an initially random placement of PS20 molecules at different concentrations. The current work investigates differences between

MA

PS20M and PS20D ester fractions and their micellar aggregates. The simulations also probe the internal reorganization of the PS20 molecule from the initial extended

TE

D

conformation, highlighting the conformational reorganization in solution due to the large number of rotatable bonds.

AC CE P

Finally, to investigate the biophysical behavior N-phenyl-1-naphthylamine (NPN), a hydrophobic fluorescence dye used in the quantification of polysorbates [40], we performed MD simulations of PS20 aggregation in the presence of NPN.

Methods

The three dimensional structures of PS20 molecules were built using MOE (MOE 2014.09; Molecular Operating Environment) for each PS20 molecule. The CHARMM general force field [41] as extension of CHARMM36 [42] was chosen for this work and PS20 atom typing and partial charges were assigned accordingly to the force field to all the atoms of PS20 molecule (Fig.1). The novel parameters of CHARMM36 led to a significant improvement in the force field in comparison with older version of CHARMM force field. For example, in some systems CHARMM36 is able to reproduce

7

ACCEPTED MANUSCRIPT a zero surface tension at the bilayer lipid density for lipids molecules [43]. Improvements are notable also in the Lennard-Jones parameters for the ester oxygen and dihedral

PT

function of the carbon chain. The geometry optimization was performed using

RI

GAUSSIAN09 [44]. The level of theory was at the B3LYP/6-31G(d) level using a tight criterion for the convergence. The geometry obtained from the B3LYP/6-31G(d) level for

SC

PS20M and PS20D were used in the following molecular mechanics simulations.

NU

Energy minimization and MD simulations were performed using NAMD version 2.9 [3]. A set of simulations was prepared for PS20M and PS20D separately, placing either

MA

molecule in a cubic box with TIP3P water model [45]. Hydrogens were added using Psfgen program in VMD [46] followed by a step with 3,000 energy minimization

TE

D

iterations carried out with NAMD version 2.9 [3] to prepare the structure. An isothermalisobaric NPT equilibration for 400 ps at 300 K was carried out with constant pressure

AC CE P

maintained at a target of 1 atm using the Nosé-Hoover-Langevin piston method [47,48]. Constant temperature was kept at 300 K using a Langevin damping coefficient of 1 ps−1. Nonbonded interactions were cut off beyond 12 Å with a smoothing function applied beyond 10 Å. Long-range electrostatics were calculated using the particle-mesh Ewald method [49]. A time step of 2 fs was used with bonded and nonbonded forces calculated every time step and particle-mesh Ewald forces calculated every other time step. System visualization, and analysis were carried out using VMD [46]. Simulations of a single molecule of PS20M, were conducted first, to determine the average conformation of the molecule in solution. This was followed by a simulation of PS20M molecules in a box at 2 different concentrations, 0.0062 M and 0.011 M, to probe micellar aggregation. A similar set of simulations was carried out for PS20D. In addition,

8

ACCEPTED MANUSCRIPT a simulation box 0.011 M with 30 molecules of PS20M and 10 molecules of PS20D has been designed to probe the effect of sorbitan monolaurate in comparison with the sorbitan

PT

dilaurate portion of PS20. A replica of PS20M system at high concentration was used to

RI

perform MD at 400 K to understand the temperature effect. The simulations of PS20M and PS20D started from an extended conformation in a non-aggregate state placing the

SC

molecule randomly in the simulation box. The system at high concentration includes up

NU

to 5.5 x 105 atoms, and it can be efficiently simulated with the use of a GPU computing cluster.

MA

The production run was 200 ns using NVT ensemble after the equilibration. Three runs were performed for each system for a total simulation time of 600 ns. The clustering

TE

D

analysis was performed based on the k-means algorithm. The algorithm starts with a global set of data from the simulations, and it returns a set of 8 clusters.

AC CE P

The adaptive biasing force (ABF) method [27] was used to estimate the free energy change along a reaction coordinate to investigate the interaction of PS20M-PS20M, PS20M-PS20D, and PS20D-PS20D, the reaction coordinate, ξ, was chosen as the distance separating the centers of mass of the two PS20 molecules. A biased force is estimated in the course of the simulation in such a way that when the force is applied to the system it results in the Hamiltonian having zero biased force along ξ. The MD simulation trajectory was generated from the bound state (ξ =10 Å) to the unbound state (ξ = 25 Å), which avoids both steric overlap in the bound state and sufficient distance to reach the unbound state based on intermolecular potentials. The trajectory was 30 ns long, and bins of 0.1 Å were used to collect the instantaneous force.

9

ACCEPTED MANUSCRIPT The structure of N-phenyl-1-naphthylamine (NPN) was optimized using GAUSSIAN09 [44], and the force field parameters were retrieved from the CHARMM general force

PT

field [41]. A simulation box at high concentration (0.011 M) of PS20M with 10

RI

molecules of NPN was prepared to analyze the interaction of NPN with PS20M. The simulations were set up and run with NAMD version 2.9 as mentioned above in the

SC

method section [3].

NU

Results

We have simulated two scenarios with 10 and 40 molecules of PS20, the latter having a

MA

concentration that is above the CMC value for PS20. The initial configuration for the PS20 mono-ester (laurate) molecule is shown in Fig. 1, after accurate initial energy

TE

D

minimization. Specifically, the simulations allow us to probe the hydrophobic effect arising from different concentrations in the simulation box, and how they result in

AC CE P

differences in PS20M and PS20D behavior in the formation of micelles. Also, we can analyze the root mean square deviations (RMSD) values, the radial distribution functions (RDFs) of PS20 molecules, the radius of gyration, and the solvent accessible surface area of the alkyl tail.

The MD simulations of the single molecule of PS20M shows that the molecule moves away from the initial extended state shown in Fig. 1, and it has a large conformational heterogeneity with a RMSD distribution that spans about 6 Å (Fig. 2). The average of the distribution is about 7.6 Å. The representative structure of the larger cluster, which is the one with a large number of members calculated with the k-means method, presents the alkyl group in such a conformational that minimizes the SASA away from the water solvent (Fig. 3). The SASA difference going from the alkyl tail exposed in solution to the

10

ACCEPTED MANUSCRIPT inward conformation is about 18.12 Å2. In the case of PS20D, the molecular dynamics simulations for the single molecule in solution showed a large conformational

PT

rearrangement. A representative conformation of a large cluster presents the two alkyls

RI

groups interacting to each other.

The simulations of the PS20M at low concentration (0.0062M) did not show the presence

SC

of a large aggregate, but we can observe only pairwise interaction between two PS20M

NU

molecules mediated by the alkyl chains that tend to stay away from the water molecules. On the other hand, the MD simulations for the high concentration scenario (above CMC)

MA

illustrate the aggregation process occurring in the PS20M. The micelle unit starts to form at about 50 ns, and then few other molecules join the micelle unit. The micelle unit has a

TE

D

roughly spherical shape with a dry hydrocarbon core, and is characterized by the presence of 18 molecules with a radius of gyration of 21 Å. In the micelle, the alkyl tail of PS20

AC CE P

flock to the interior of the micelle, excluding water molecules from the highly hydrophobic environment (Fig. 4). The large polar head of the PS20 resides in the outer surface of the micelle surrounded by and interacting with water molecules in the solvent. Few intermolecular hydrogen bonds stabilized the polar region of PS20. Overall, the hydrophobic regions of the PS20 minimize the SASA, lying in the middle of the micelle unit where van der Waals contacts also contribute to the alkyl chain interactions. As general mechanism we can observe the two PS20 molecules interact with each other creating an initial hydrophobic environment where other PS20 molecules diffusing in solution can eventually approach to rearrange their conformations, placing their hydrophobic tails in close proximity with those of the initial two molecules. Once the

11

ACCEPTED MANUSCRIPT micelle is formed, it is quite stable with the PS20 in the micelle rarely exchanging with those in aqueous solution for the duration of the ~200 ns simulation.

PT

The simulations of the PS20D show a much stronger hydrophobic effect than PS20M;

RI

indeed we can observe the formation of a micelle unit after 20 ns instead of 50 ns for PS20M. This result is consistent across the three replicate simulations, which utilized

SC

random initial placement of the PS20 molecules in the simulation box. The micelle unit

NU

involves 25 PS20D molecules, and shows a high degree of compactness in comparison with the micelle of PS20M (Fig. 4). The hydrocarbon tails of the two arms rearrange in a

MA

conformation with the tails pointing towards the core of spherical micelle. The density of these tails in the core is higher than in the micelle unit of PS20M. The radius of gyration

TE

D

of the micellar aggregate for both PS20D and PS20M are fairly similar, about 22 Å, but the PS20D core is denser. The length of the hydrocarbon chain impacts the diameter of

AC CE P

the micelle (steric hindrance of the head region is a factor as well), but the molecular packing is due to the presence of two tails instead of one. The simulations including both PS20M and PS20D molecules in the box also show that a mixed micelle can be formed with a radius of gyration of about 23 Å with 21 molecules. As shown in Fig. 5 the RDF at 300 K of the PS20 for all atoms presents a main peak at 1.1 Å, a second peak is located at 1.5 Å, and a third peak above 100 RDF value is present at 2.2 Å. The first peak remarks the short-range interaction due to Van der Waals forces of the hydrocarbon group present in the interior of the micelle. It is a typical first peak appearing in the molecular simulations of Lennard-Jones fluids. The second main peak shows a lower RDF value but still above 100 at 1.5 Å. The third main peak is located 2.2 Å. Essentially no structure can be seen beyond 4 Å. Instead, in the simulation at 400 K

12

ACCEPTED MANUSCRIPT (useful for determining qualitative trends) we can note how the main peak at 1.1 Å shows a higher RDF value. This suggests an aggregation process that could be driven by an

PT

entropic effect (e.g. increase in water entropy). A less pronounced increase can be

RI

noticed for the second and the third main peaks as well. Similarly, the RDF plot for the 400 K simulations shows no structure beyond 4 Å. Overall, the occurrence of the three

SC

main peaks highlights that micelle exhibits a well-defined structure.

NU

We can measure the total number of water-water hydrogen bonds in the simulations at low and high concentration (Table 1). The normalized number per molecule of water-

MA

water hydrogen bonds for the low concentration of PS20M is 0.6105. Instead, the high concentration PS20 simulation shows 0.6154. This increase may seem slight, however

TE

D

could be quite significant if one considers the change in the amount of water interacting with the surfactant molecules. This agrees with previous work pointing out that the total

AC CE P

number of hydrogen bonds in the system always increase in the process of formation of micelles [50]. This trend reflects significant differences in the hydrogen bond network due the high number of water molecules included in the simulation box in case of the formation of the micelle.

Aside from water-water hydrogen bonds, the hydrogen bonds patterns also give a contribution to stabilize energetically the hydrophilic head of PS20 as in Fig. 6. Time evolution of the distance ξ during the 30 ns of ABF simulation of the set of simulations of PS20M-PS20M, PS20D-PS20M, and PS20D-PS20D is achieved. In the first 10 ns, the two PS20 molecules progressively get away from each other as the energy barriers implied by the breaking of van der Waals interactions are sampled and overcome.

13

ACCEPTED MANUSCRIPT The molecules then get close to each other, and ξ begins to fluctuate freely, providing evidences that the free energy barrier is efficiently overcome during the simulations.

PT

As shown in Fig. 7 the free energy profile for the dissociation of PS20M-PS20M presents

RI

a barrier of about 19 kcal/mol. The profile for PS20M-PS20D has a higher barrier (about 21 kcal/mol) than PS20M-PS20M. A larger barrier in comparison with the other two sets

NU

energy are located at ξ = 10.1 Å in all three cases.

SC

has to be overcome for the dissociation of PS20D-PS20D (~28 kcal/mol). Minima of free

Lastly, the MD simulations of NPN in the simulation box with PS20M and PS20D

MA

confirm that this small molecule interacts with PS20. In particular, after the formation of the micelle we can find the majority of NPN molecules (13 out of 15) located in the core

TE

D

of the micelle unit. The NPN molecules are packed together with the hydrocarbon atoms in the hydrophobic region of the micelle (Fig. 8). The aryl-alkyl interactions stabilize the

AC CE P

NPN molecule in the hydrophobic core, and these are weak interactions dominated by London dispersion forces. A close-up view as in Fig. 9 shows how the NPN molecules are located in the middle of the micelle core, and the water molecules are away from this core. As shown in Fig. 10 the distance of water molecules from the core of the micelle unit is typically greater than 8 Å for PS20M and 10 Å for PS20D. The peak of this distribution is at 25 Å for PS20M and 28 Å for PS20D, and these distributions confirm that the water molecules mostly reside away from the dry core of the micelle, and in the case of PS20D the waters are shifted away from the core of the micelle.

Discussion The conformational heterogeneity of PS20M and PS20D reflects the presence of many conformational states in solution, and confirms the complexity of the dynamics that

14

ACCEPTED MANUSCRIPT involved many degrees of freedom for PS20. Importantly, the conformations of the PS20 in solution migrate away from the initial extended state, and assume a compact

PT

conformation according the minima of the energy potential. The dynamics of the

RI

molecule affects the formation of the micelle as well. In particular, the aggregation behavior of PS20 at high concentration agrees well with the biophysical trends observed

SC

experimentally in which the micelle unit forms at and above the CMC value [51]. The

NU

free surfactant concentration fluctuated substantially throughout the dynamic equilibrium simulation, and N=3 simulation replicates were not statistically sufficient to determine an

MA

accurate CMC value, unfortunately. However, the aggregate number found in the micellar aggregates were much more stable across over the duration of the simulation and

TE

D

its replicate runs. Also at sufficiently low concentrations, the simulations show no formation of micelles. The radius of gyration of the micelle unit of about 21-22 Å agrees

AC CE P

with experimental findings coming from scattering techniques, which is about 24 Å [2]. Since we have used only the laurate as fatty acid fraction of the PS20 tail, which roughly represents the ~40-60 % of the possible mixture of fatty acids tails, differences found in scattering experiments are expected. However, the MD simulations provided insightful trends and interesting details on the structure of the micelle at atomic resolution. The hydrocarbon atoms oriented toward the core of the micelle confirms the extent of the hydrophobic effect for the PS20 molecules. The entropy drives the formation of the micelle due to the fact that the solvation free energy for the micelle unit is lower than the sum of solvation free energies for individual particles. Although away from the micelle core, the structure and the dynamics of water molecules around the micelle play a crucial role. Mixed micelles made of PS20M and PS20D randomly placed in the simulation box

15

ACCEPTED MANUSCRIPT are possible, leading to the conclusion that in solution, a single micelle can be composed by different fractions of PS20.

PT

The hydrogen bond network of water changes from low concentration, in which no

RI

micelle is formed, to the high PS20 concentration case. The number of hydrogen bonds increases when micelle units form. The large size of the micelle makes it impossible to

SC

keep the same hydrogen bonding patterns that can be accommodated for small

NU

hydrophobic molecules. These results also highlight the crucial role of waters in the PS20 aggregation as we can see the hydrogen bond network significantly change when the

MA

micelle is formed increasing the number of hydrogen bonds per water molecule. This evidence has been reported also using temperature-jump Monte Carlo simulations

TE

D

directed to analyze the formation of hydrogen bonds and micelle dynamics [52]. Our results are in line with experimental and calculations previously performed to tackle

AC CE P

the hydrophobic effect problem [30,31,35,53]. The RDF plots show also how the high temperature increases the main peak of the RDF which corresponds to an increase of attractive forces for PS20M and is in line with experimental evidences that CMC decreases with increased temperature. This suggests that the micellar aggregation could be partly driven by an entropic effect (e.g. related to water rearrangement). Overall, the hydrophobic forces involved in the formation of the micelle show a much greater magnitude. Structural evaluation of the system revealed that penetration of water molecules into the micelle was restricted to the head group region, leaving a relatively water-free hydrocarbon core. When the system is at low temperature, there is a reduction in the entropic term of the free energy, which leads to a lower enthalpy through an increase in

16

ACCEPTED MANUSCRIPT the total number of hydrogen bonds in the system involving the hydrophilic head portion. The ABF calculations provide evidence for a stronger free energy of reversible

PT

association for PS20D than PS20M, which agrees with other results showing that in the

RI

simulation of PS20D at high concentration the micelle unit forms in a shorter timescale (20 ns) than for PS20M (50 ns). Furthermore, this can help explain why in general the

SC

CMC decreases with the length of the hydrocarbon chain of the tail group, given the fact

NU

that the favorable free energy component of the solvation free energy derives mostly from the transfer of the hydrophobic portion from the polar solvent environment to the

MA

core of the micelle unit [54].

The NPN molecule simulations shed insight on the behavior of this particular

TE

D

hydrophobic molecule that has been used in the biochemical assay to estimate the CMC value of micelle formation for PS20. Our results highlight the strong interaction of the

AC CE P

NPN with PS20 molecules, in particular the NPN molecules after the formation of the micelle are located in the hydrophobic core of the micelle binding with aryl-alkyl interactions to the hydrocarbon atoms and staying away from the water molecules in the dry core. The strong interaction NPN-PS20 can explain the 6-fold difference in the CMC values using the surface tension measurement compared to the NPN assay as shown in previous experiments [51].

Conclusions The MD simulations of the PS20 molecule show a large conformational change in solution starting from an extended configuration, due to the multiple possible conformational states. The most probable of a PS20 molecule on its own tends to minimize exposure of its hydrophobic tail. The micellar aggregation of PS20 molecules

17

ACCEPTED MANUSCRIPT in solution is governed by the hydrophobic effect assuming a particular conformational arrangement that minimizes the solvation free energy. The length of the hydrocarbon tail

PT

of PS20, the number of its tails, and the temperature rise can affect the extent of the

RI

hydrophobic effect and the density of the micelle unit. The micelle structures found in the MD simulations done here are quite stable, with practically all hydrophobic tail atoms

SC

located at the center of the micelle. The hydrophilic heads of the micelle can form

NU

hydrogen bonds among each other and constitute a conformational steric and entropic constraint on the size of the micelle unit. Different variants of PS20 (PS20M, and PS20D)

MA

can form mixed micelle units. The NPN dye interacts with the hydrophobic tail of PS20, and during the formation of the micelle is found in the hydrophobic core of the micelle

TE

D

forming aryl-alkyl type of interaction. Water molecules tend to stay away from the dry core of the micelle unit, and less water is found in PS20D micelles than in PS20M

AC CE P

micelles. The computational work presented here advances the understanding of the thermodynamic forces that drive the biophysical behavior of PS20 in solution, providing atomistic details of micelle structure and NPN interaction with the micelle.

18

ACCEPTED MANUSCRIPT

References

AC CE P

TE

D

MA

NU

SC

RI

PT

[1] B.A. Kerwin, Polysorbates 20 and 80 used in the formulation of protein biotherapeutics: Structure and degradation pathways, J. Pharm. Sci. 97 (2008) 2924-2935. [2] C.C. Ruiz, J.A. Molina-Bolivar, J. Aguiar, G. MacIsaac, S. Moroze, R. Palepu, Effect of ethylene glycol on the thermodynamic and micellar properties of Tween 20, Colloid. Polym. Sci. 281 (2003) 531-541. [3] J.C. Phillips, R. Braun, W. Wang, J. Gumbart, E. Tajkhorshid, E. Villa, C. Chipot, R.D. Skeel, L. Kale, K. Schulten, Scalable molecular dynamics with NAMD, J. Comput. Chem. 26 (2005) 1781-1802. [4] D.P. Tieleman, D. van der Spoel, H.J.C. Berendsen, Molecular dynamics simulations of dodecylphosphocholine micelles at three different aggregate sizes: Micellar structure and chain relaxation, J. Phys. Chem. B 104 (2000) 6380-6388. [5] S. Bogusz, R.M. Venable, R.W. Pastor, Molecular dynamics simulations of octyl glucoside micelles: Structural properties, J. Phys. Chem. B 104 (2000) 54625470. [6] L. Cunha-Silva, J.J.C. Teixeira-Dias, Aqueous solution inclusion of the nonionic surfactant C12E4 in beta-cyclodextrin: Implications of micellization in stoichiometry determination and model calculations, J. Inclusion Phenom. Macrocycl. Chem. 43 (2002) 127-131. [7] S.J. Marrink, D.P. Tieleman, A.E. Mark, Molecular dynamics simulation of the kinetics of spontaneous micelle formation, J. Phys. Chem. B 104 (2000) 12165-12173. [8] T. Lazaridis, B. Mallik, Y. Chen, Implicit solvent simulations of DPC micelle formation, J. Phys. Chem. B 109 (2005) 15098-15106. [9] R.A. Karjiban, M. Basri, M.B.A. Rahman, A. Salleh, Structural Properties of Nonionic Tween80 Micelle in Water Elucidated by Molecular Dynamics Simulation, 2nd International Conference on Chemistry and Chemical Process (Icccp 2012) 3 (2012) 287-297. [10] C.D. Bruce, M.L. Berkowitz, L. Perera, M.D.E. Forbes, Molecular dynamics simulation of sodium dodecyl sulfate micelle in water: Micellar structural characteristics and counterion distribution, J. Phys. Chem. B 106 (2002) 3788-3793. [11] S. Storrn, S. Jakobtorweihen, I. Smirnova, A.Z. Panagiotopoulos, Molecular Dynamics Simulation of SDS and CTAB Micellization and Prediction of Partition Equilibria with COSMOmic, Langmuir 29 (2013) 11582-11592. [12] X. Tang, P.H. Koenig, R.G. Larson, Molecular Dynamics Simulations of Sodium Dodecyl Sulfate Micelles in Water—The Effect of the Force Field, J. Phys. Chem. B 118 (2014) 3864-3880. [13] A.F. de Moura, L.C.G. Freitas, Molecular dynamics simulation of the sodium octanoate micelle in aqueous solution, Chem. Phys. Lett. 411 (2005) 474-478.

19

ACCEPTED MANUSCRIPT

AC CE P

TE

D

MA

NU

SC

RI

PT

[14] A.F. de Moura, L.C.G. Freitas, Molecular dynamics simulation of the sodium octanoate micelle in aqueous solution: Comparison of force field parameters and molecular topology effects on the micellar structure., Braz. J. Phys. 34 (2004) 64-72. [15] E. Gallicchio, M. Lapelosa, R.M. Levy, The Binding Energy Distribution Analysis Method (BEDAM) for the Estimation of Protein-Ligand Binding Affinities., J. Chem. Theory Comput. 6 (2010) 2961-2977. [16] M. Lapelosa, C.F. Abrams, A computational study of water and CO migration sites and channels inside myoglobin, J. Chem. Theory Comput. 9 (2013) 12651271. [17] M. Lapelosa, E. Gallicchio, R.M. Levy, Conformational Transitions and Convergence of Absolute Binding Free Energy Calculations, J. Chem. Theory Comput. 8 (2012) 47-60. [18] M.S. Head, J.A. Given, M.K. Gilson, “Mining Minima”:  Direct Computation of Conformational Free Energy, J. Phys. Chem. A 101 (1997) 1609-1618. [19] L. Maragliano, A. Fischer, E. Vanden-Eijnden, G. Ciccotti, String method in collective variables: Minimum free energy paths and isocommittor surfaces, J. Chem. Phys. 125 (2006). [20] W.L. Jorgensen, C. Ravimohan, Monte-Carlo Simulation of Differences in FreeEnergies of Hydration, J. Chem. Phys. 83 (1985) 3050-3054. [21] A. Laio, M. Parrinello, Escaping free-energy minima, Proc. Natl. Acad. Sci. U. S. A. 99 (2002) 12562-12566. [22] T. Huber, A.E. Torda, W.F. Vangunsteren, Local Elevation - a Method for Improving the Searching Properties of Molecular-Dynamics Simulation, J. Comput. Aided Mol. Des. 8 (1994) 695-708. [23] D. Hamelberg, J.A. McCammon, Standard free energy of releasing a localized water molecule from the binding pockets of proteins: Double-decoupling method, J. Am. Chem. Soc. 126 (2004) 7683-7689. [24] D.L. Mobley, J.D. Chodera, K.A. Dill, On the use of orientational restraints and symmetry corrections in alchemical free energy calculations, J. Chem. Phys. 125 (2006). [25] P.A. Kollman, I. Massova, C. Reyes, B. Kuhn, S.H. Huo, L. Chong, M. Lee, T. Lee, Y. Duan, W. Wang, O. Donini, P. Cieplak, J. Srinivasan, D.A. Case, T.E. Cheatham, Calculating structures and free energies of complex molecules: Combining molecular mechanics and continuum models, Acc. Chem. Res. 33 (2000) 889897. [26] M. Lapelosa, T.W. Patapoff, I.E. Zarraga, Modeling of protein–anion exchange resin interaction for the human growth hormone charge variants, Biophys. Chem. 207 (2015) 1-6. [27] E. Darve, A. Pohorille, Calculating free energies using average force, J. Chem. Phys. 115 (2001) 9169-9183. [28] J. Hénin, A. Pohorille, C. Chipot, Insights into the Recognition and Association of Transmembrane α-Helices. The Free Energy of α-Helix Dimerization in Glycophorin A, J. Am. Chem. Soc. 127 (2005) 8478-8484.

20

ACCEPTED MANUSCRIPT

AC CE P

TE

D

MA

NU

SC

RI

PT

[29] J. Comer, J.C. Gumbart, J. Hénin, T. Lelièvre, A. Pohorille, C. Chipot, The Adaptive Biasing Force Method: Everything You Always Wanted To Know but Were Afraid To Ask, J. Phys. Chem. B 119 (2015) 1129-1151. [30] D. Chandler, Interfaces and the driving force of hydrophobic assembly, Nature 437 (2005) 640-647. [31] L.R. Pratt, A. Pohorille, Hydrophobic Effects and Modeling of Biophysical Aqueous Solution Interfaces, Chem. Rev. 102 (2002) 2671-2692. [32] N. Giovambattista, C.F. Lopez, P.J. Rossky, P.G. Debenedetti, Hydrophobicity of protein surfaces: Separating geometry from chemistry, Proc. Natl. Acad. Sci. U. S. A. 105 (2008) 2274-2279. [33] J.A. Morrone, J. Li, B.J. Berne, Interplay between Hydrodynamics and the Free Energy Surface in the Assembly of Nanoscale Hydrophobes, J. Phys. Chem. B 116 (2012) 378-389. [34] J. Kyte, The basis of the hydrophobic effect, Biophys. Chem. 100 (2002) 193203. [35] L.F. Scatena, M.G. Brown, G.L. Richmond, Water at hydrophobic surfaces: Weak hydrogen bonding and strong orientation effects, Science 292 (2001) 908912. [36] T.W. Randolph, L.S. Jones, Surfactant-protein interactions, Pharm. Biotechnol. 13 (2002) 159-175. [37] V.B. Fainerman, R. Miller, E.V. Aksenenko, A.V. Makievski, J. Krägel, G. Loglio, L. Liggieri, Effect of surfactant interfacial orientation/aggregation on adsorption dynamics, Adv. Colloid Interface Sci. 86 (2000) 83-101. [38] R. Braun, D.M. Engelman, K. Schulten, Molecular Dynamics Simulations of Micelle Formation around Dimeric Glycophorin A Transmembrane Helices, Biophys. J. 87 (2004) 754-763. [39] J. Huang, A.D. MacKerell, CHARMM36 all-atom additive protein force field: Validation based on comparison to NMR data, J. Comput. Chem. 34 (2013) 2135-2145. [40] M. Khossravi, Y.H. Kao, R.J. Mrsny, T.D. Sweeney, Analysis methods of polysorbate 20: A new method to assess the stability of polysorbate 20 and established methods that may overlook degraded polysorbate 20, Pharm. Res. 19 (2002) 634-639. [41] K. Vanommeslaeghe, E. Hatcher, C. Acharya, S. Kundu, S. Zhong, J. Shim, E. Darian, O. Guvench, P. Lopes, I. Vorobyov, A.D. MacKerell, CHARMM General Force Field: A Force Field for Drug-Like Molecules Compatible with the CHARMM All-Atom Additive Biological Force Fields, J. Comput. Chem. 31 (2010) 671-690. [42] R.B. Best, X. Zhu, J. Shim, P.E. Lopes, J. Mittal, M. Feig, A.D. Mackerell, Jr., Optimization of the additive CHARMM all-atom protein force field targeting improved sampling of the backbone phi, psi and side-chain chi(1) and chi(2) dihedral angles, J. Chem. Theory Comput. 8 (2012) 3257-3273. [43] R.W. Pastor, A.D. Mackerell, Jr., Development of the CHARMM Force Field for Lipids, J. Phys. Chem. Lett. 2 (2011) 1526-1532. [44] M.J. Frisch, G.W. Trucks, H.B. Schlegel, G.E. Scuseria, M.A. Robb, J.R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G.A. Petersson, H. Nakatsuji, M. Caricato, 21

ACCEPTED MANUSCRIPT

AC CE P

TE

D

MA

NU

SC

RI

PT

X. Li, H.P. Hratchian, A.F. Izmaylov, J. Bloino, G. Zheng, J.L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J.A. Montgomery, J.E. Peralta, F. Ogliaro, M. Bearpark, J.J. Heyd, E. Brothers, K.N. Kudin, V.N. Staroverov, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J.C. Burant, S.S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J.M. Millam, M. Klene, J.E. Knox, J.B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R.E. Stratmann, O. Yazyev, A.J. Austin, R. Cammi, C. Pomelli, J.W. Ochterski, R.L. Martin, K. Morokuma, V.G. Zakrzewski, G.A. Voth, P. Salvador, J.J. Dannenberg, S. Dapprich, A.D. Daniels, Farkas, J.B. Foresman, J.V. Ortiz, J. Cioslowski, D.J. Fox, Gaussian 09, Revision B.01, Wallingford CT, 2009. [45] W.L. Jorgensen, J. Chandrasekhar, J.D. Madura, R.W. Impey, M.L. Klein, Comparison of Simple Potential Functions for Simulating Liquid Water, J. Chem. Phys. 79 (1983) 926-935. [46] W. Humphrey, A. Dalke, K. Schulten, VMD: Visual molecular dynamics, J. Mol. Graphics 14 (1996) 33-38. [47] G.J. Martyna, D.J. Tobias, M.L. Klein, Constant-Pressure Molecular-Dynamics Algorithms, J. Chem. Phys. 101 (1994) 4177-4189. [48] S.E. Feller, Y.H. Zhang, R.W. Pastor, B.R. Brooks, Constant-Pressure MolecularDynamics Simulation - the Langevin Piston Method, J. Chem. Phys. 103 (1995) 4613-4621. [49] T. Darden, D. York, L. Pedersen, Particle Mesh Ewald - an N.Log(N) Method for Ewald Sums in Large Systems, J. Chem. Phys. 98 (1993) 10089-10092. [50] G. Heinzelmann, W. Figueiredo, M. Girardi, Micellar dynamics and water-water hydrogen-bonding from temperature-jump Monte Carlo simulations, Chem. Phys. Lett. 550 (2012) 83-87. [51] A. Patist, S.S. Bhagwat, K.W. Penfield, P. Aikens, D.O. Shah, On the measurement of critical micelle concentrations of pure and technical-grade nonionic surfactants, J. Surfactants Deterg. 3 (2000) 53-58. [52] G. Heinzelmann, W. Figueiredo, M. Girardi, Micellar dynamics and water–water hydrogen-bonding from temperature-jump Monte Carlo simulations, Chem. Phys. Lett. 550 (2012) 83-87. [53] E. Gallicchio, M.M. Kubo, R.M. Levy, Enthalpy-entropy and cavity decomposition of alkane hydration free energies: Numerical results and implications for theories of hydrophobic solvation, J. Phys. Chem. B 104 (2000) 6271-6285. [54] L. Maibaum, A.R. Dinner, D. Chandler, Micelle formation and the hydrophobic effect, J. Phys. Chem. B 108 (2004) 6778-6781.

22

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

MA

Fig. 1: The three dimensional structure of PS20M with five repeats of ethylene oxide subunits for each arm for a total of 20 (left). The partial charges of the sorbitan head

AC CE P

TE

D

group used in the molecular simulations (right).

23

AC CE P

TE

D

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

Fig. 2: Time evolution (top) and distribution (bottom) of all atoms RMSD for PS20M with respect to the initial conformation for the molecular simulations of the single molecule in solution.

24

TE

D

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

Fig. 3: Representative snapshot of the larger cluster of a single PS20M obtained with k-

AC CE P

means clustering.

25

SC

RI

PT

ACCEPTED MANUSCRIPT

NU

Fig. 4: Representative structures of the micelle for PS20M (left) and PS20D (right). The structure is represented in stick with the oxygen atoms in red and the carbon in sandy

AC CE P

TE

D

MA

brown. The hydrophobic tails are colored in green.

26

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

AC CE P

TE

D

MA

Fig. 5: RDF of PS20M atoms from the high concentration simulation at 300 K and 400 K.

27

D

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

AC CE P

TE

Fig. 6: Hydrogen bonds between two hydrophilic heads in the micelle unit of PS20M.

28

AC CE P

TE

D

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

29

ACCEPTED MANUSCRIPT

Fig. 7: Free energy profiles for PS20M-PS20M (top), PS20M-PS20D (middle), and

PT

PS20D-PS20D (bottom) reversible association. Beyond point 20 Å we can consider it not interacting. In order to estimate the variation three ABF simulations of 5 ns were carried

AC CE P

TE

D

MA

NU

SC

RI

out independently for each PMF, starting from the current PMF curves.

30

D

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

TE

Fig. 8: Structural representation of NPN molecules interacting with the micelle unit. The NPN molecules are interacting in the core of the micelle. The NPN molecules are

AC CE P

represented in green, and the PS20 atoms are in sandy brown with stick representation.

31

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

Fig. 9: Close up view of the NPN molecule (green) interacting with hydrocarbon atoms in

D

the middle of the micelle core with water molecules with the oxygen depicted in red and

AC CE P

TE

the hydrogen colored in white located away from the center of the micelle unit.

32

AC CE P

TE

D

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

Fig. 10: Distributions of distance of the waters molecule from the center of mass of the micelle in PS20M and PS20D.

33

ACCEPTED MANUSCRIPT Table 1: Comparison of the hydrogen bonds (HB) network in the PS20M simulation at high concentration with the low concentration simulation. The total average is calculated

PT

averaging out all the frames of the MD trajectory. The normalized HB number is

Total average number of HB

Normalized HB number

5.0893e+04 1.0426e+05

0.6105 0.6154

SC

Concentration

AC CE P

TE

D

MA

NU

Low Conc. High Conc.

RI

calculated per water molecule.

34

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

AC CE P

TE

D

MA

Graphical abstract

35

ACCEPTED MANUSCRIPT

AC CE P

TE

D

MA

NU

SC

RI

PT

Highlights  Molecular simulations show micellar aggregation of PS20  Molecular simulations probe the hydrophobic effect  Simulations provide atomistic details of the micelle unit for PS20M and PS20D  Differences in the hydrogen bond network are present for PS20M and PS20D

36