Optimizing Protein–Polymer Interactions in a Poly(ethylene glycol

Aug 14, 2018 - Optimizing Protein–Polymer Interactions in a Poly(ethylene glycol) Coarse-Grained ... when these two groups of materials integrate to...
0 downloads 0 Views 4MB Size
Download PDF
Article Cite This: J. Phys. Chem. B XXXX, XXX, XXX−XXX

pubs.acs.org/JPCB

Optimizing Protein−Polymer Interactions in a Poly(ethylene glycol) Coarse-Grained Model Farhad Ramezanghorbani,† Ping Lin,† and Coray M. Colina*,†,‡ †

Department of Chemistry and ‡Department of Material Science and Engineering, University of Florida, Gainesville, Florida 32611, United States

Downloaded via KAROLINSKA INST on August 19, 2018 at 02:50:48 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

S Supporting Information *

ABSTRACT: Increasing demand for hybrid materials that merge the synthetic and biological areas in drug industries requires in-depth knowledge of the individual components and their contributions to these complexes. Coarse-grained (CG) models developed for proteins and polymers exist, yet there is a lack of understanding of the cross interactions when these two groups of materials integrate to build a complex. In this work, we characterized the nonbonded interactions between poly(ethylene glycol) (PEG) and amino acids in a Martini CG model utilizing state-of-the-art quantum mechanics calculations of interaction energies. The parameter set proposed, was validated by assessing the polymer density in the vicinity of individual amino acids obtained from available all-atomistic molecular dynamic simulations of plasma proteins. Our results revealed the necessity of protein−polymer interaction parameterization at the CG level to avoid overestimation of polymer association when employing other PEG models within the Martini framework.

■

INTRODUCTION Protein-based drugs have been extensively used in pharmaceutical industries for their enhanced stability and controllable activity.1,2 PEGylation, a process of covalent attachment of a substrate to biocompatible poly(ethylene glycol) (PEG), an amphiphilic polymer, soluble both in polar and apolar solvents, is extensively employed for a broad range of applications in biomedical industries.1−3 Protein interactions with amphiphilic polymers enhance the solubility of the complex, and prevent antibody recognition, while expanding the hydrodynamic volume and therefore increasing the drug’s circulating lifetime.4 The kidney clearance molecular weight (MW) cutoff is estimated to be 70 kDa,2 therefore, for most of the plasma proteins, with MW around 65 kDa, the circulating lifetime of a protein−polymer conjugate could increase significantly by employing low MW PEG (∼5 kDa). This could reduce one of the problems in using high MW PEG in drug complexes, the degradation of PEG in the human body.3 The dependencies of the retention time, solubility, and effective volume on the MW of conjugates have been investigated using gel permeation chromatography and electrophoresis techniques.5,6 However, rationalizing the effect of the underlying protein−polymer interactions, in the exposure of proteins functional epitopes and enhanced hydrodynamic behavior of the complex, are more difficult to access by experimental methods. Molecular dynamics (MD) simulations are common in silico methods extensively utilized to study these hybrid compounds.7,8 For example, Settanni et al.9,10 © XXXX American Chemical Society

characterized the interactions between PEG and several plasma proteins using atomistic MD simulations. PEG−protein interactions were compared for individual amino acid type in terms of preferential binding coefficients and the binding ratio of cosolvent (PEG) to solvent (water). The amino acid specific binding affinity was shown to hold within the error bar across several plasma proteins. PEG was shown to have more favorable interactions with apolar residues than polar and charged ones.9,10 Consequently, PEG is not totally inert, and its interaction with proteins depends on the amino acid composition of the surface. Although atomistic simulations provide a realistic description of molecular systems, their application is limited by the size and the timescale of the phenomena under study. The flexibility of PEG chains as well as the high molecular mass of protein−polymer conjugates introduce a lack of sampling in high-resolution atomistic simulations due to the large number of particles and slow dynamics of the system. Alternatively, lower-resolution coarse-grained (CG) models provide enhanced sampling of conformational phase space while maintaining an accurate thermodynamic and structural view of many proteins11,12 and polymers.13−15 These models are developed by mapping a group of atoms, residues, or part of the system into a single interaction site (also called bead or Received: June 4, 2018 Revised: July 24, 2018

A

DOI: 10.1021/acs.jpcb.8b05359 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 1. Martini CG nonpolarizable mapping for individual amino acids, labeled with the corresponding Martini types. The backbone bead type for GLY, ALA, PRO, and other amino acids depends on the location of the amino acid (i.e., secondary structure’s dependency). Different types assigned for (1) coil, bend, or free, (2) helix, (3) helix C/N terminus, and (4) β-strand or turn configurations are as follows: GLY: (P5, N0, Nd/Na, Nda), ALA: (P4, C5, N0/N0, N0), PRO: (P4, C5, N0/Na, N0), and for the rest of the backbone beads: (P5, N0, Nd/Na, Nda).

properties calculated from higher-resolution atomistic simulations.14,15 These PEG CG models were often employed in CG studies together with other compounds, such as proteins, lipids, etc.7,22,23 However, these PEG CG models were not specifically parameterized to reproduce the nonbonded interactions with other proteins and often overestimate or underestimate the association of PEG chains with them.18 The Martini forcefield is limited to only 10 discrete pair interaction strengths (epsilon of Lennard-Jones potential) and 18 different bead definitions (subcategories Q, P, N, and C) to describe complex atomic groups.12 These inherent approximations of the Martini forcefield and in general, coarse-graining, demand further changes when other molecules coexist with PEG polymers in the system under study. Failing to correctly address the cross interactions may result in inaccurate prediction of structures, activity, and stability of protein− polymer systems. For years, PEG interactions with proteins and peptides were often approximated by taking into account the excluded volume effect. For example, Hamed et al.22 investigated the effect of conjugation site on the stability and assembly mechanism of PEGylated four helix bundle using the CG model developed by Lee et al.14 The nonbonded interactions between PEG and the protein in their work were approximated by an entropic repulsion, assuming PEG was inert, and thus only resembling the excluded volume effect. However, recent studies have shown that PEG interactions with protein are complex.9,10 The nonbonded interaction parameters in the PEG CG model developed by Lee et al.14 were adjusted to reproduce

super atom). This process allows for longer integration timesteps and enhances the dynamics by reducing the number of particles, the degrees of freedom, and consequently smoothing out the free-energy landscape. The Martini forcefield is one of the well-established CG forcefields employed for a wide range of chemical compounds.12−21 A combination of bottom-up and top-down approaches was employed, by utilizing atomistic simulations and experimental partitioning free energies, to provide a hierarchical parametrization of the so-called “Martini building blocks”. Typically, the bonded interaction parameters are derived from high-resolution atomistic simulations or crystallographic databases, whereas experimental free energies of transfer between water and organic solvents are employed for tuning the nonbonded interactions. Every 2, 3, or 4 heavy atoms, with their associated hydrogens, are mapped into one CG bead. However, due to the impractical choice of 4-to-1 mapping to describe ringlike structures and small units, “Stype” CG beads with 2-to-1 or 3-to-1 mapping were introduced by scaling down the original effective radii and interaction strengths. Four main categories are proposed for CG bead types in Martini: Q, for charged, P for polar, N for intermediate, and C for apolar compounds, each with subcategories that distinguish the level of polarity or hydrogen-bonding capacity (donor−acceptor). The PEG CG models developed within the Martini framework exhibit reasonable agreement with experimental densities, neutron scattering data, and MW dependency of structure as well as the structural and thermodynamic B

DOI: 10.1021/acs.jpcb.8b05359 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

secondary and tertiary structures of proteins in this forcefield due to the lack of hydrogen bond directionality. Thus, the secondary and tertiary structures are enforced by an elastic network of harmonic bonds. The elastic bond force constants were set to 1000 kJ mol−1 nm−2 for intrahelical and 250 kJ mol−1 nm−2 for interhelical distances with lower and upper elastic bond cutoffs of 0.5 and 0.9 nm, respectively. The PEG CG model is based on a three-to-one mapping of a C−O−C repeat unit to an S-type Martini bead, as displayed in Figure 2. In this model, one bead type is used for simplicity.

the radius of gyration, obtained by atomistic simulations in explicit solvent, and the experimental densities of low MW PEG. Bonded degrees of freedom in this model were parameterized using probability distributions obtained from atomistic simulations performed with the CHARMM ether forcefield.24 The bonded degrees of freedom in this PEG CG model allow for stable CG MD simulations with no more than 8 fs integration timestep, which is significantly shorter than the timestep typically used for proteins with the Martini forcefield (20−25 fs). This model will be referred to as sN0 throughout this work. Although the model reproduces the structural and thermodynamic properties of PEG in both atomistic simulations and experiments, it fails in reproducing the thickness of the PEG layer on PEGylated lipids, compelling the authors to refine the nonbonded interaction parameters.18 In another study, Rossi et al.15 developed a PEG CG model, in which the bonded interactions were parameterized to reproduce the probability distribution of bonded degrees of freedom obtained from atomistic simulations with the GROMOS forcefield.25,26 The nonbonded interactions were assigned in a process of matching the CG octanol−water partitioning free energies of building blocks to the experimental values. However, the absence of a set of parameters in the “Martini matrix” to satisfy these free energies compelled the authors to introduce an intermediate level among predefined interaction types in the matrix. We refer to this model as sP0 in this study. In this work, we evaluate the applicability of the protein− polymer nonbonded interaction parameters in the sP0 and sN0 models and propose a new bottom-up approach to characterize these interactions in a PEG CG model using state-of-the-art quantum mechanics (QM) energy calculations. These interactions were further validated by reproducing the local density of PEG in the vicinity of various plasma proteins obtained by atomistic simulations. The explicit solvent CG MD simulations of various plasma proteins and free short PEG polymers revealed that the previous PEG models14,15 (sN0, sP0) present significant deviations from the atomistic simulation results in terms of PEG’s coordination and affinity toward individual amino acids. The overestimation of these properties is resolved with the introduction of a set of nonbonded interactions for the new PEG model (sNP). Importantly, the optimized PEG CG model (sNP) is not inert and it shows different affinities with different amino acids.

Figure 2. PEG coarse-grained representation for four PEG repeat units. Three heavy atoms (C−O−C) were mapped into a CG bead highlighted in blue.

The interaction strength (Lennard-Jones epsilon value) between S-type pairs is reduced to 75% of their original value, and the effective interaction size between them (sigma value) is 0.43 nm instead of the original 0.47 nm. The bonded degrees of freedom (bonds, angles, and dihedrals) together with nonbonded polymer−polymer and polymer−water interactions were adopted from the model developed by Rossi et al.15 This model allows stable CG MD simulations with integration timestep of 20−25 fs, in line with other Martini models. The nonbonded interactions between PEG beads and amino acid CG beads were evaluated based on three different models: (i) the standard Martini sN0 bead type assigned to the PEG bead as suggested by Lee et al.;14 (ii) the sP0 bead type assigned to the PEG bead developed by Rossi et al.;15 and (iii) assignment of a new sNP bead type to PEG model, as described in detail below. The interactions between the new sNP bead type and the existing protein CG beads were estimated using QM calculations. The chemical characteristics (polarity, composition, etc.) of the PEG beads and the association energies between PEG beads and models of amino acids building blocks were thus taken into account effectively. QM Estimation of Nonbonded Interactions between CG Model Beads. To assign the appropriate parameters for the nonbonded interactions between PEG CG beads and protein CG beads, we used QM calculations to determine the relative strengths between CG bead pairs. Dimethoxyethane (DME) can be seen as a poly(ethylene oxide) or PEG dimer and is often used as the smallest model for them. Therefore, we used the interaction energies obtained by QM as a guide to assign the nonbonded interactions between PEG CG beads and protein CG beads. The calculations were performed using Gaussian 09 D.01.27 Geometry optimizations and frequency calculations were carried out at the MP2/6-31G(d) level of theory using the SMD solvation model to represent the aqueous environment.28 The selection of the reference structures for CG beads is provided in Table 1. The interaction energies (ΔE) between model molecules, A and B, were obtained using eq 1.

■

METHODS Coarse-Grained Models for Proteins and PEG. The CG models for bovine serum albumin (BSA), human serum albumin (HSA), and apo-human serum transferrin (apo-hTF) plasma proteins were generated by using the Martinize python script16 with the PDB ids 4F5S, 1AO6, and 2HAV, respectively. In these crystal structures, BSA and apo-hTF do not present missing residues. HSA does have four missing residues in the N terminal and three missing residues in the C terminal that were not reconstructed. The CG representation of individual amino acids and the Martini bead type assignment are shown in Figure 1. It is important to mention that alanine, glycine, proline, and the backbone beads can adopt different Martini types depending on the secondary structure environment in the protein’s crystal structure. The bonded and nonbonded interactions of the protein CG models are described in detail in the Martini forcefield v2.2.12,16,17 Additional bias is required to preserve the

ΔE = E(A + B) − E(A) − E(B) C

(1) DOI: 10.1021/acs.jpcb.8b05359 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

that enhances kinetics, the effective time sampled with Martini CG models is estimated12 to be typically 4-fold longer than the mentioned trajectory length (2 μs). The equilibrated box sizes together with the total number of particles (PEG, water, protein, and ions) for simulation of BSA, HSA, and apo-hTF are provided in Table S1. A new interaction set, NP, is introduced to describe the new PEG bead type. The mapping definition for sNP is equivalent to the sP0 and sN0 models of Lee et al.14 and Rossi et al.15 In addition to the bonded interaction parameters, the set of parameters describing the nonbonded interactions between sNP−sNP and sNP−P4 (water) was taken from the sP0 model. All parameters for the new sNP model can be downloaded from http://forcefield-database.org. The nonbonded interactions between our polymer model, sNP, and the other 18 Martini bead types were assigned and validated by calculating a quantity relevant to excess of PEG monomers in the local environment of the protein during CG MD trajectories. This quantity measures the stoichiometric number ratio of PEG monomers and water beads within a certain cutoff distance from each amino acid type, where the stoichiometric number can be obtained using eq 2. Normalizing this ratio with its corresponding bulk value (eq 3) allows transferability of calculated ratios among different simulation conditions. For consistency with calculations obtained from atomistic simulations,9,10 we refer to this ratio as PEGWAA (eq 4).

Table 1. Reference Structures for QM Calculations of CG Bead Interaction Energies with DME CG bead

model molecules (SMILES format)

Qda Qd Qa Q0 P5 P4 P3 P2 P1 Nda Nd Na N0 C5 C4 C3 C2 C1

[NH3+][CH2][CH2][OH] [NH3+][CH2][CH2][CH3] [CH3][CH2]C(O)[O−] [CH3][NH+]([CH3])[CH3] [NH2]C(O)[CH3] [OH][CH2][CH2][OH] [CH3][NH][CH]O [CH3][CH2][OH] [CH3][CH2][CH2][OH] [CH3][CH2][CH2][CH2][OH] [CH3][CH2][CH2][NH2] [CH3]C(O)[CH3] [CH3][CH2]O[CH3] [CH3][CH2][CH2][SH] [CH2][CH][CH][CH2] [CH3][CH][CH][CH3] [CH3][CH2][CH3] [CH3][CH2][CH2][CH3]

corresponding amino acid groups Lys side chain Asp, Glu side chain Asn, Gln side chain, backbone

Ser, Thr side chain

Cys side chain

Val side chain Ile side chain

Since each DME contains two PEG CG beads, the calculated interaction energies should be proportional to the interaction between a single PEG CG bead and each model CG bead listed in Table 1. Although this procedure is not an accurate assessment, such interaction energies provide useful information to guide the ranking and assignment of the nonbonded interactions between CG beads. CG MD Simulation Details. Coarse-grained MD simulations were performed with the GROMACS 2016.4 package.29 The temperature was set to 300 K using the velocity rescaling stochastic algorithm30 with a coupling time constant of 1.0 ps. The pressure was kept at 1.0 bar using the Parrinello−Rahman barostat31,32 with a coupling time constant of 12.0 ps and a compressibility factor of 3 × 10−4 bar−1. A timestep of 20 fs with a leapfrog MD integrator was used, whereas the neighbor list was updated every 400 fs using a Verlet neighbor list scheme33 to impose a reduced straight cutoff of 1.1 nm for nonbonded interaction calculations. In the Verlet scheme, the neighbor list is automatically built based on the energy drift tolerance. In this work, a Verlet buffer tolerance of 0.005 kJ mol−1 ps−1 was selected. Long-range electrostatics were calculated with the reaction-field algorithm34 using a cutoff of 1.1 nm and a relative dielectric constant of 15. Each protein CG structure was energy-minimized prior to solvation with nonpolarizable Martini water beads, each bead equivalent to four water molecules. Then, 64 short polymer chains, each with 4 monomers (PEG4), were randomly inserted into a cubic box. To achieve neutral systems, 17, 15, and 2 sodium ion beads were respectively added to the boxes containing BSA, HSA, and apo-hTF by random substitution of water beads. A secondary relaxation was performed in the presence of solvent, followed by a 5 ns isothermal−isobaric (NPT) equilibration while restraining bead positions in each protein. Initial structures for production MD were obtained by performing further NPT equilibration for 10 ns while the additional restraints were released. Production runs were performed in the NPT ensemble for 2 μs, with a 20 fs integration timestep. However, it is important to indicate that due to the relatively smooth underlying free-energy landscape

Nj =

∫0

Rc

g j,AA (r )4πr 2 dr

ij NPEG yz ij N yz jj zz = lim R c →∞jjj PEG zzz jj N zz j Nwater z k water { bulk k { PEGWAA =

NPEG/Nwater (NPEG/Nwater)bulk

(2)

(3)

(4)

Where Nj is the coordination number of particle j in the vicinity of each amino acid (AA) type, defined as distances closer than Rc = 0.7 nm, the cutoff, from the AA. Integration of the pair correlation function for particles i and j (gij) result in a coordination number with a certain cutoff. The first minima along the radial distribution function of polymer−amino acids in the system were chosen for the cutoff in our calculations. The larger the binding affinity between amino acid and PEG, the higher the PEGWAA ratio. The average and standard deviation for PEGWAA is calculated using 4 × 500 ns blocks of MD trajectory. Graphics were generated using the visual molecular dynamics (VMD) package.35

■

RESULTS AND DISCUSSION QM-Guided CG Parameterization for Amino Acid− PEG Interactions. Methoxyethane, the PEG monomer, is slightly more polar than the N0 Martini CG bead type used to describe it. This is reflected in the reproduction of experimental partitioning free energies between apolar−polar phases in the Martini CG forcefield, in which the N0 bead type slightly underestimates the experimental values of the reference compound, methoxyethane.12 Furthermore, the necessity of optimizing dipole−dipole, and dipole-induced dipole interactions between PEG and other bead types, has been already demonstrated by the enhancement of the PEG−water and PEG−PEG interactions in the original sN0 model14 in a PEGylated lipid study.18 Alternatively, even though our D

DOI: 10.1021/acs.jpcb.8b05359 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B analysis reveals that the sP0 model15 overestimates the polar− polar interaction of PEG with amino acids, the model maintains similar interaction parameters to the sN0 model for polymer−water and polymer−polymer interactions. These effective nonbonded interaction strengths are εself = 0.75 × 4.5 kJ mol−1 (Martini level II) for both models and εsP0−water = 4.25 kJ mol−1, exceeding the corresponding sN0−water value by 0.25 kJ mol−1. Although both models reasonably reproduce PEG’s structural and thermodynamic properties, we implemented the bonded interactions and nonbonded polymer− polymer, polymer−water interaction parameters as developed in the sP0 model.15 As mentioned before, the bonded interaction parameters in the sP0 model enable the use of longer timesteps (20−25 fs) compatible with other Martini molecules, in contrast with the 8 fs timestep of the sN0 model. QM calculations of interaction energies were used to assist in the assignment of a new bead type description for PEG monomers to specifically include targeted PEG−amino acid nonbonded interactions. This new bead type is an intermediary level between N0 and P0, and will be referred to as NP throughout this work. The assignment of pairwise nonbonded interaction levels, the most rigorous task in the parametrization of any CG model, become a tractable task by using the Martini discrete interaction matrix (Table S2). Here, we limited the nonbonded interaction assignment for NP−X pairs to supraattractive O, attractive I, almost attractive II, semiattractive III, intermediate IV, almost intermediate V, and semirepulsive VI Martini levels. In Figure 3 are shown the QM interaction energies calculated for the 18 reference functional groups (correspond-

ing to each Martini bead type) and a DME molecule (a PEG dimer). The configuration dependency (trans/cis) of the dipole moment25,36 in PEG is taken into account by employing QM calculations on DME instead of the PEG monomer. DME can adopt a cis configuration in the vicinity of nonpolar groups and a trans configuration while interacting with polar functional groups to respectively benefit from hydrophobic effects and dipole−dipole interactions. Amino acid interactions cannot directly be determined from PEGWAA ratios or QM energies. Many body effects are not taken into account for the PEGWAA ratio calculations, i.e., the interactions do not scale linearly to the rates or their exponentials, and the contributions from neighboring amino acids also contribute to their values. Additionally, the inherent loss of entropy in CG modeling needs to be compensated by enthalpic protein−polymer interactions that are respectively weaker than QM energies. In short, the QM energy values in Figure 3 were employed indirectly to assign relative interaction strengths in the C, N, P, and Q Martini bead type categories. The lower limit for interaction of PEG with other Martini types was set to semirepulsive VI (ε = 2.7 kJ mol−1) as in all of the intermediate N-type beads interacting with C1 (rightmost column in Table S2). The sNP interactions with apolar bead types were assigned to the same semirepulsive VI level, and were gradually shifted to level IV for C5, consistent with the decrease observed in the calculated QM energies. Since PEG can adopt dihedral distributions that minimize or maximize the dipole moment, the polymer can have preferential interactions toward both apolar and polar compounds. The hydrogen bond acceptor nature of DME is reflected in the QM energy calculations; lower energies are associated with hydrogen bond donors, and thus intermediate bead type N, with donor (Nd), or donor-and-acceptor (Nda) hydrogen bond capabilities were treated with a level III parameter. Interaction with P4 (Martini nonpolarizable water) was fixed to the numerical average of values for level II and III, to match the corresponding value in the sP0 model.15 Polar beads P3 and P5 were also set to the same value of P4 based on the QM calculations. In the charged region of the Martini matrix, Q, the interaction levels were assigned according to the order of QM energies, with the highest level, O for Q0 and the lowest, IV for Qa. Our final parameter set for the nonbonded interactions (NP) together with N0 and P0 parameter levels

Figure 3. QM interaction energies calculated using a set of 18 reference compounds corresponding to each Martini bead type.

Figure 4. CG interaction levels for N0 and P0 models compared with NP. Interaction levels of a Martini water bead (P4) with all 18 Martini types are plotted in blue. The order of these values when interacting with different bead types and the difference with respect to the interaction with water reveals the competition of cosolvent−solvent in the system. E

DOI: 10.1021/acs.jpcb.8b05359 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 5. PEG−amino acid contact distribution, normalized by number of trajectory frames for sNP, sN0, and sP0 models interacting with BSA. Disagreement for affinity patches, and their location across the protein surface may cause inaccurate prediction of structural and thermodynamic properties of protein−polymer systems in the cases of sN0 and sP0.

fraction of solvent accessible surface residues that are exposed to PEG monomers during the MD trajectory (Table S1). The percentage coverage of ∼100% in all cases ensures all of the amino acids on the protein surface were explored by the polymer CG model. The nonbonded parameters in our model were compared with the sN0 and sP0 PEG CG models in the reproduction of spots on the protein surface that have high preferential interactions with PEG. The relative strength of these highaffinity patches and their location on the protein surface were monitored using the average density of PEG monomers on the surface of BSA. In Figure 5 are shown these patches on the BSA’s surface, obtained from averages over 2 μs CG MD trajectories. On the left side of Figure 5, the distribution of the number of contacts that occurred within 0.7 nm of a polymer− amino acid distance is shown. Variation in the peak locations and the height of the distributions indicate that high-affinity patches on the protein surface vary as a function of the CG model used. The optimized sNP model has a relatively lower binding affinity toward BSA when compared against the sN0 and sP0 models. This inconsistency between CG models could consequently affect the prediction of protein activity by incorrect masking of functional groups. Additionally, overestimating the PEG’s association to protein surfaces may result in underestimated hydrodynamic radius, retention time, and

are shown together with the original Martini matrix in Table S2. The resulting cosolvent−solvent competition is also demonstrated by a comparison of the effective interaction strengths given in Figure 4. The difference in effective interaction of each individual Martini type X with water, P4−X, and different levels used for the PEG parameterization, NP−X, N0−X, and P0−X, indicates that in nearly all of the models, apolar interactions are favored. The discrepancy arises in the extent of this difference for C-type, and the order of it in the N, P, and Q regions of the parameter set. In contrast with other models, our parametrization process distinguishes between charged types by taking into account the difference in energy of these building blocks interacting with DME. Contacts between Bovine Serum Albumin and PEG4. CG solutions of short PEG chains and different plasma proteins (BSA, HSA, and apo-hTF) were prepared to reproduce the corresponding PEGWAA ratios provided by the atomistic simulations of Settanni et al.9,10 The reference allatomistic PEGWAA values were obtained from Figure 3 of reference 10. In all cases, 64 short PEG chains of length 4 monomers (PEG4) in both atomistic and CG MD simulations diffuse faster than long PEG chains, allowing for a better sampling of the protein surface. The sampling quality was monitored by percentage coverage (C%), defined as the F

DOI: 10.1021/acs.jpcb.8b05359 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

secondary and tertiary structure have been found to be intact upon conjugation and adsorption of different MW PEG onto their surface. The validation of the sNP model is thus further extended to other plasma proteins HSA and apo-hTF. We compare the CG residue-specific affinities for our PEG model with the atomistic simulation data provided in the work of Settanni et al.9,10 for HSA, BSA, and apo-hTF to validate the set of parameters used for describing PEG’s interaction with different amino acids. Figure 7 illustrates the correlation between PEGWAA ratios

circulating lifetime of protein−polymer conjugates that function as drugs. PEG−Amino Acid Interactions in BSA. The parameters for the sNP model were validated against all-atomistic PEGWAA ratios previously calculated by Settanni et al.9,10 for various plasma proteins. Figure 6 provides a comparison

Figure 6. PEGWAA ratios of CG models compared to atomistic simulations of Settanni et al.9,10 values for each amino acid type in BSA. The local density of the sNP model reproduces the results provided by the atomistic simulations, whereas the sN0 and sP0 models overestimate the ratios. The bars for CG ratios indicate the standard deviation of 4 block averages over 2 μs simulations. The dashed line corresponds to PEGWAA = 1, where the amino acid affinities toward PEG and water are similar. For the values above the dashed line, PEG is more favorable.

between our model sNP and the sN0 and SP0 models and with respect to the atomistic PEGWAA ratios for amino acids in BSA. A PEGWAA ratio greater than one indicates the higher affinity of amino acid toward PEG with respect to water. The deviation of this quantity in Figure 6 is associated with the surface accessibility of the amino acids. Exposed charged and hydrophilic amino acids cover more than 70% of the protein surface (Table S3), thus they have smaller deviations from PEGWAA mean value; the deviation is larger for the apolar residues (e.g., M) due to their rare appearance on the protein surface. Aside from the chemical nature of the amino acids, their availability also affects the precision of the PEGWAA ratio calculated for them. For example, BSA with 583 amino acids, has only two Tryptophan (W) residues; this results in higher standard deviation calculated for its PEGWAA ratio. Our coarse-grained model reproduces the atomistic PEGWAA ratio, whereas the other two CG models (sN0 and sP0) are overestimating PEG’s binding to the protein surface. Moreover, our model successfully reproduced the all-atomistic9,10 high-density regions of PEG near the protein surface (Figure S1). Plasma Protein−PEG4 Interactions in PEGWAA. BSA and HSA are capable of transporting drugs and fatty acids while maintaining the osmotic pressure of the blood.37−39 apohTF is used as a target specific drug delivery agent to brain tissues through blood.40,41 Therefore, the effect of PEG− plasma protein interactions in modulating drug delivery and reducing their clearance rate from the body makes these proteins important candidates for protein−polymer interaction validation in this work. PEGylation is shown to have no significant effect on the structure of these plasma proteins in multiple experimental37,39,42 studies in which protein’s

Figure 7. CG PEGWAA correlation with ratios calculated with respect to the atomistic MD simulations for plasma proteins: BSA, HSA, and apo-hTF. The corresponding error bars for atomistic and CG calculation are included only for ALA. The dashed line represents perfect agreement with the atomistic simulations. Solid lines represent linear fit to the simulation data. The linear correlation of our sNP model (black) indicates the accuracy of the parameters used to describe the PEG−amino acid interactions. The relatively lower rates in apo-hTF emphasize the protein-dependent nature of the PEG− amino acid interactions.

calculated in atomistic and CG MD simulations for BSA, HSA, and apo-hTF. The accuracy of the interaction parameters for the sNP model is highlighted by the proximity of the CG PEGWAA results with the reference line for each amino acid type (perfect agreement with the atomistic simulations of Settanni et al.9,10 is indicated by the dashed line). The reduced overall PEG−protein affinity of the new CG model (sNP) compared to sN0 and sP0 models is mainly attributed to the difference in the interactions with C* bead G

DOI: 10.1021/acs.jpcb.8b05359 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

from the new PEG CG model are in good agreement with the results from atomistic simulations for BSA, HSA, and apohTF plasma proteins, the similarity of PEGWAA across number of plasma proteins9,10 suggests that this model can be used for other members in this group of proteins.

types. The percentage of solvent accessible surface area (SASA) for different bead types existing in the protein CG models (Figure 1) is provided in Table S3. More than 70% of the surface in BSA, HSA, and apo-hTF is covered with polar (P) and charged bead types (Q), however, the overall affinity of PEG in our model is mainly affected by the remaining ∼30% apolar (C) and intermediate (N) bead types. Although, Qaand Qd-charged groups are the dominant surface-exposed Martini beads, they predominantly favor water (Figure 4). Alternatively, C* bead types interact favorably with PEG models (Figure 4). Although fewer C* beads exist on the protein surface, their weaker interactions (with sNP model) affect the overall protein−PEG affinity, leading to significantly improved PEGWAA ratio as shown in Figure 7. Although the PEGWAA ratio was shown above to be consistent across several plasma proteins,9,10 the amino acid composition of the protein surface could affect the calculated values by the explicit incorporation of cross interactions and many body effects. However, the results here presented confirm the transferability of our parameters for HSA, BSA, and apo-hTF. The PEGWAA ratio for individual amino acid type differs for the different proteins, since the surface composition of the protein dictates the binding properties of PEG toward individual amino acids. As a result, apo-hTF disfavors PEG interactions relative to HSA and BSA. With 75.8% sequence identity to HSA, BSA has a different surface amino acid composition,38 which also explains the difference in their PEGWAA ratios.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.8b05359. Information related to the size of simulations and surface residue coverage obtained; Martini nonbonded interaction matrix including set of parameters for PEG CG bead type (including the N0, P0, and new NP types); percentage SASA of individual bead types on protein surface for BSA, HSA, and apo-hTF; PEG’s density map on BSA surface; PEGWAA ratios of CG models compared to atomistic references (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]fl.edu. Phone: +1 (352) 294 3488. ORCID

Farhad Ramezanghorbani: 0000-0002-7545-4416 Ping Lin: 0000-0002-6141-5424 Coray M. Colina: 0000-0003-2367-1352

■

Notes

CONCLUSIONS We presented a new set of nonbonded potential parameters that describe PEG−protein interactions in a CG model within the Martini framework. The bonded degrees of freedom and polymer−polymer, polymer−water interactions were adopted from the model developed by Rossi et al.15 The strength of the nonbonded Lennard-Jones interactions between the new PEG model (sNP) and the remaining Martini bead types was adjusted by comparing QM energies between DME, a PEG dimer, and the corresponding reference structures of the Martini CG bead types. Our model was validated by reproducing specific affinities between amino acids and PEG, from the high-resolution atomistic simulations of Settanni et al.9,10 The local density ratio of the PEG CG monomers to water, in the vicinity of amino acids, was found in agreement with the ratios obtained from atomistic simulations for various plasma proteins BSA, HSA, and apo-hTF. We found that in contrast with our model, other PEG CG models developed within the Martini framework (Lee et al.,14 and Rossi et al.15) overestimate the protein−polymer interactions and present different amino acid specificities. This is due to the limited levels of interactions in the Martini matrix that could represent the correct balance between polymer−polymer, polymer− water, and polymer−protein interactions. Consequently, employing the previous PEG CG models (Lee et al.,14 and Rossi et al.15) in complex with protein CG models may provide less accurate predictions of the properties, such as protein’s activity, hydrodynamic radius, etc. Our nonbonded potential parameter set for PEG (sNP) is reliable for protein−polymer applications within the Martini framework. Moreover, our bottom-up approach is not limited to the models developed within the Martini forcefield and could be applied to other CG models with proven representability. Although, the PEGWAA ratios calculated

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors acknowledge the financial support from University of Florida Preeminence Initiative and the support for computational resources provided by UF Research Computing and Information Technology Center at University of Florida.

■

REFERENCES

(1) Alconcel, S. N. S.; Baas, A. S.; Maynard, H. D. FDA-approved poly(ethylene glycol)−protein conjugate drugs. Polym. Chem. 2011, 2, 1442. (2) Caliceti, P.; Veronese, F. M. Pharmacokinetic and biodistribution properties of poly(ethylene glycol)-protein conjugates. Adv. Drug Delivery Rev. 2003, 55, 1261−1277. (3) Veronese, F. M.; Pasut, G. PEGylation, successful approach to drug delivery. Drug Discovery Today 2005, 10, 1451−1458. (4) Shu, J. Y.; Panganiban, B.; Xu, T. Peptide-polymer conjugates: from fundamental science to application. Annu. Rev. Phys. Chem. 2013, 64, 631−657. (5) Sherman, M. R.; Williams, L. D.; Saifer, M. G. P.; French, J. A.; Kwak, L. W.; Oppenheim, J. J. In Poly(ethylene glycol) Chemistry and Biological Applications; Harris, J. M., Zalipsky, S., Eds.; American Chemical Society: Washington, DC, 1997; Vol. 680, pp 155−169. (6) Nakaoka, R.; Tabata, Y.; Yamaoka, T.; Ikada, Y. Prolongation of the serum half-life period of superoxide dismutase by poly(ethylene glycol) modification. J. Controlled Release 1997, 46, 253−261. (7) Lee, H. Molecular modeling of PEGylated peptides, dendrimers, and single-walled carbon nanotubes for biomedical applications. Polymers 2014, 6, 776−798. (8) Stanzione, F.; Jayaraman, A. Hybrid atomistic and coarse-grained molecular dynamics simulations of polyethylene glycol (PEG) in explicit water. J. Phys. Chem. B 2016, 120, 4160−4173. (9) Settanni, G.; Zhou, J.; Suo, T.; Schöttler, S.; Landfester, K.; Schmid, F.; Mailänder, V. Protein corona composition of poly(ethylene glycol)- and poly(phosphoester)-coated nanoparticles H

DOI: 10.1021/acs.jpcb.8b05359 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B correlates strongly with the amino acid composition of the protein surface. Nanoscale 2017, 9, 2138−2144. (10) Settanni, G.; Zhou, J.; Schmid, F. Interactions between proteins and poly(ethylene-glycol) investigated using molecular dynamics simulations. J. Phys.: Conf. Ser. 2017, 921, No. 012002. (11) Ramezanghorbani, F.; Dalgicdir, C.; Sayar, M. A multi-state coarse grained modeling approach for an intrinsically disordered peptide. J. Chem. Phys. 2017, 147, No. 094103. (12) Marrink, S. J.; Risselada, H. J.; Yefimov, S.; Tieleman, D. P.; de Vries, A. H. The MARTINI force field: coarse grained model for biomolecular simulations. J. Phys. Chem. B 2007, 111, 7812−7824. (13) Panizon, E.; Bochicchio, D.; Monticelli, L.; Rossi, G. MARTINI coarse-grained models of polyethylene and polypropylene. J. Phys. Chem. B 2015, 119, 8209−8216. (14) Lee, H.; de Vries, A. H.; Marrink, S.-J.; Pastor, R. W. A coarsegrained model for polyethylene oxide and polyethylene glycol: conformation and hydrodynamics. J. Phys. Chem. B 2009, 113, 13186−13194. (15) Rossi, G.; Fuchs, P. F. J.; Barnoud, J.; Monticelli, L. A coarsegrained MARTINI model of polyethylene glycol and of polyoxyethylene alkyl ether surfactants. J. Phys. Chem. B 2012, 116, 14353− 14362. (16) de Jong, D. H.; Singh, G.; Bennett, W. F. D.; Arnarez, C.; Wassenaar, T. A.; Schäfer, L. V.; Periole, X.; Tieleman, D. P.; Marrink, S. J. Improved parameters for the MARTINI coarse-grained protein force field. J. Chem. Theory Comput. 2013, 9, 687−697. (17) Monticelli, L.; Kandasamy, S. K.; Periole, X.; Larson, R. G.; Tieleman, D. P.; Marrink, S.-J. The MARTINI coarse-grained force field: Extension to proteins. J. Chem. Theory Comput. 2008, 4, 819− 834. (18) Lee, H.; Pastor, R. W. Coarse-grained model for PEGylated lipids: effect of PEGylation on the size and shape of self-assembled structures. J. Phys. Chem. B 2011, 115, 7830−7837. (19) Menichetti, R.; Kanekal, K. H.; Kremer, K.; Bereau, T. In silico screening of drug-membrane thermodynamics reveals linear relations between bulk partitioning and the potential of mean force. J. Chem. Phys. 2017, 147, No. 125101. (20) Frederix, P. W. J. M.; Patmanidis, I.; Marrink, S. J. Molecular simulations of self-assembling bio-inspired supramolecular systems and their connection to experiments. Chem. Soc. Rev. 2018, 47, 3470− 3489. (21) Shen, H.; Moustafa, I. M.; Cameron, C. E.; Colina, C. M. Exploring the dynamics of four RNA-dependent RNA polymerases by a coarse-grained model. J. Phys. Chem. B 2012, 116, 14515−14524. (22) Hamed, E.; Ma, D.; Keten, S. Effect of polymer conjugation site on stability and self-assembly of coiled coils. BioNanoScience 2015, 5, 140−149. (23) Hamed, E.; Xu, T.; Keten, S. Poly(ethylene glycol) conjugation stabilizes the secondary structure of α-helices by reducing peptide solvent accessible surface area. Biomacromolecules 2013, 14, 4053− 4060. (24) Vorobyov, I.; Anisimov, V. M.; Greene, S.; Venable, R. M.; Moser, A.; Pastor, R. W.; MacKerell, A. D. Additive and classical drude polarizable force fields for linear and cyclic ethers. J. Chem. Theory Comput. 2007, 3, 1120−1133. (25) Fuchs, P. F. J.; Hansen, H. S.; Hünenberger, P. H.; Horta, B. A. C. A GROMOS parameter set for vicinal diether functions: properties of polyethyleneoxide and polyethyleneglycol. J. Chem. Theory Comput. 2012, 8, 3943−3963. (26) Horta, B. A. C.; Fuchs, P. F. J.; van Gunsteren, W. F.; Hü nenberger, P. H. New interaction parameters for oxygen compounds in the GROMOS force field: Improved pure-liquid and solvation properties for alcohols, ethers, aldehydes, ketones, carboxylic acids, and esters. J. Chem. Theory Comput. 2011, 7, 1016−1031. (27) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; et al. Gaussian 09, revision D.01; Gaussian, Inc.: Wallingford, CT, 2009.

(28) Marenich, A. V.; Cramer, C. J.; Truhlar, D. G. Universal solvation model based on solute electron density and on a continuum model of the solvent defined by the bulk dielectric constant and atomic surface tensions. J. Phys. Chem. B 2009, 113, 6378−6396. (29) Abraham, M. J.; Murtola, T.; Schulz, R.; Páll, S.; Smith, J. C.; Hess, B.; Lindahl, E. GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers. SoftwareX 2015, 1−2, 19−25. (30) Bussi, G.; Donadio, D.; Parrinello, M. Canonical sampling through velocity rescaling. J. Chem. Phys. 2007, 126, No. 014101. (31) Parrinello, M.; Rahman, A. Polymorphic transitions in single crystals: A new molecular dynamics method. J. Appl. Phys. 1981, 52, 7182. (32) Parrinello, M.; Rahman, A. Crystal structure and pair potentials: A molecular-dynamics study. Phys. Rev. Lett. 1980, 45, 1196−1199. (33) Páll, S.; Hess, B. A flexible algorithm for calculating pair interactions on SIMD architectures. Comput. Phys. Commun. 2013, 184, 2641−2650. (34) Onsager, L. Electric moments of molecules in liquids. J. Am. Chem. Soc. 1936, 58, 1486−1493. (35) Humphrey, W.; Dalke, A.; Schulten, K. VMD: visual molecular dynamics. J. Mol. Graphics 1996, 14, 33−38. (36) Smith, G. D.; Borodin, O.; Bedrov, D. A revised quantum chemistry-based potential for poly(ethylene oxide) and its oligomers in aqueous solution. J. Comput. Chem. 2002, 23, 1480−1488. (37) Plesner, B.; Fee, C. J.; Westh, P.; Nielsen, A. D. Effects of PEG size on structure, function and stability of PEGylated BSA. Eur. J. Pharm. Biopharm. 2011, 79, 399−405. (38) Bujacz, A. Structures of bovine, equine and leporine serum albumin. Acta Crystallogr., Sect. D: Biol. Crystallogr. 2012, 68, 1278− 1289. (39) Akbarzadehlaleh, P.; Mirzaei, M.; Mashahdi-Keshtiban, M.; Shamsasenjan, K.; Heydari, H. Pegylated human serum albumin: review of pegylation, purification and characterization methods. Adv. Pharm. Bull. 2016, 6, 309−317. (40) Jain, A.; Chasoo, G.; Singh, S. K.; Saxena, A. K.; Jain, S. K. Transferrin-appended PEGylated nanoparticles for temozolomide delivery to brain: in vitro characterisation. J. Microencapsulation 2011, 28, 21−28. (41) Kim, T. H.; Jo, Y. G.; Jiang, H. H.; Lim, S. M.; Youn, Y. S.; Lee, S.; Chen, X.; Byun, Y.; Lee, K. C. PEG-transferrin conjugated TRAIL (TNF-related apoptosis-inducing ligand) for therapeutic tumor targeting. J. Controlled Release 2012, 162, 422−428. (42) Ferebee, R.; Hakem, I. F.; Koch, A.; Chen, M.; Wu, Y.; Loh, D.; Wilson, D. C.; Poole, J. L.; Walker, J. P.; Fytas, G.; et al. Light scattering analysis of mono- and multi-PEGylated bovine serum albumin in solution: Role of composition on structure and interactions. J. Phys. Chem. B 2016, 120, 4591−4599.

I

DOI: 10.1021/acs.jpcb.8b05359 J. Phys. Chem. B XXXX, XXX, XXX−XXX