Effects of Polymer Conjugation on Hybridization Thermodynamics of

Article pubs.acs.org/JPCB

Effects of Polymer Conjugation on Hybridization Thermodynamics of Oligonucleic Acids Ahmadreza F. Ghobadi† and Arthi Jayaraman*,†,‡ †

Department of Chemical and Biomolecular Engineering, Colburn Laboratory, University of Delaware, 150 Academy Street, Newark, Delaware 19711, United States ‡ Department of Material Science and Engineering, University of Delaware, Newark, Delaware 19711, United States S Supporting Information *

ABSTRACT: In this work, we perform coarse-grained (CG) and atomistic simulations to study the effects of polymer conjugation on hybridization/melting thermodynamics of oligonucleic acids (ONAs). We present coarse-grained Langevin molecular dynamics simulations (CG-NVT) to assess the effects of the polymer flexibility, length, and architecture on hybridization/ melting of ONAs with different ONA duplex sequences, backbone chemistry, and duplex concentration. In these CGNVT simulations, we use our recently developed CG model of ONAs in implicit solvent, and treat the conjugated polymer as a CG chain with purely repulsive Weeks−Chandler−Andersen interactions with all other species in the system. We find that 8− 100-mer linear polymer conjugation destabilizes 8-mer ONA duplexes with weaker Watson−Crick hydrogen bonding (WC Hbonding) interactions at low duplex concentrations, while the same polymer conjugation has an insignificant impact on 8-mer ONA duplexes with stronger WC H-bonding. To ensure the configurational space is sampled properly in the CG-NVT simulations, we also perform CG well-tempered metadynamics simulations (CG-NVT-MetaD) and analyze the free energy landscape of ONA hybridization for a select few systems. We demonstrate that CG-NVT-MetaD simulation results are consistent with the CG-NVT simulations for the studied systems. To examine the limitations of coarse-graining in capturing ONA− polymer interactions, we perform atomistic parallel tempering metadynamics simulations at well-tempered ensemble (AAMetaD) for a 4-mer DNA in explicit water with and without conjugation to 8-mer poly(ethylene glycol) (PEG). AA-MetaD simulations also show that, for a short DNA duplex at T = 300 K, a condition where the DNA duplex is unstable, conjugation with PEG further destabilizes DNA duplex. We conclude with a comparison of results from these three different types of simulations and discuss their limitations and strengths.

I. INTRODUCTION Oligonucleic acid (ONA) hybridization is a thermoreversible process where two single strands of ONAs form a duplex via interstrand hydrogen bonding of complementary bases and stacking interactions between neighboring bases on the same strand.1 Upon increasing the system temperature, the duplex melts into two single-stranded ONAs. Due to the reversibility and specificity of the hybridization, ONA, and in particular DNA, is the basis of many applications in gene sequencing,2,3 molecular recognition and biosensing,4−6 DNA switches,7 selfassembly of nanoscale materials,8−11 and DNA origami.12−14 To control and program these diverse applications, much work has been done to understand how various ONA design parameters, such as duplex length, sequence, concentration, chemistry, G···C content, base mismatches, salt concentration, and hydrophilic and hydrophobic surface chemistries, impact ONA hybridization/melting.15−21 Additionally, the quest for cheaper analogues of the universally used DNA has led to alternative synthetic chemistries that also exhibit thermorever© XXXX American Chemical Society

sibility and specificity with increased chemical diversity and reduced biological stability of DNA.22−30 To understand how these new chemistries within ONA duplexes impact their melting/hybridization, in a recent work,31 we studied the effects of backbone flexibility, backbone charge, and nucleobase distance on duplex hybridization thermodynamics. We showed using molecular simulations that decreasing ONA backbone flexibility, ONA backbone charge density, and intrastrand nucleobase distance increases ONA duplex stability. In many applications where DNA or new alternative ONAs are used, such as biosensing,32,33 gene delivery,34−38 ONA based networks/gels for tissue engineering,39−45 and ONA based nanostructures,46−51 ONAs are conjugated to other biocompatible polymers, like poly(ethylene glycol) (PEG). Therefore, there is a critical need to understand how the Received: July 12, 2016 Revised: August 26, 2016

A

DOI: 10.1021/acs.jpcb.6b06970 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 1. Schematic of our coarse-grained model for a representative polymer-conjugated ONA duplex with sequence CGCGCGCG. Naturally occurring bases C and G are shown in red and blue, respectively. All backbone (BB) beads are of the same type but colored differently for visual depiction of the sequence. Hydrogen bond (HB) sites are different for each natural base but are all colored in yellow. Polymer (PL) beads do not have a HB site, and are colored in gray.

thermodynamics of generic ONAs in implicit solvent and explicit ions.31 In this section, we present only the key features of our ONA CG model, and direct the reader to the original publication for further details. In our CG model, each nucleotide is represented with one coarse-grained backbone (BB) bead (Figure 1). Nucleobase distance and single strand flexibility are controlled by introducing harmonic BB−BB bonds and BB−BB−BB angles, respectively. Depending on the ONA backbone chemistry, the BB beads is negatively charged as in DNA or neutral as in peptide nucleic acid73 (PNA) or click-nucleic acid45 (CNA). The interstrand Watson−Crick hydrogen bonds (WC H-bonds) between complementary bases are captured by introducing one electrostatically neutral hydrogen bonding (HB) site on each BB bead. There are four types of HB sites (A, T, C, and G) to mimic the four naturally occurring bases. To include stacking interactions between neighboring nucleobases on the same strand, a harmonic dihedral potential is defined for every HB−BB− BB−HB 1−4 interaction.31 In this work, we introduce a new electrostatically neutral bead type, PL, to construct the polymer chain for the polymer−ONA conjugate (Figure 1). Harmonic bond potentials are used to connect adjacent PL−PL and PL− BB beads, where the equilibrium bond distance and the harmonic potential constant are similar to that of the BB−BB bonds.31 Harmonic angle potentials are defined for PL−PL−PL angles to control the polymer flexibility. The BB−BB−PL and BB−PL−PL angles have ONA-like (kangle = 0−30 ε/rad2) and polymer-like (kangle = 0−30 ε/rad2) flexibilities, respectively. No other dihedral angle potential is introduced in the polymer chain. The last bead type is IN to represent explicit monovalent ions, with charge valency of +1 or −1 for positive and negative ions, respectively. The units of length, energy, and mass are defined as σ = 0.6 nm, ε = 0.1 kcal/mol, and m = 42.0 g/mol, respectively. The diameters of BB, HB, PL, and IN beads are 1.0σ, 0.3σ, 1.0σ, and 0.7σ, respectively. The masses of BB, HB, PL, and IN beads are 3.0m, 1.0m, 1.0m, and 1.0m, respectively. The size and mass of BB beads roughly represent the size and mass of a single DNA nucleobase. Because our CG model is a generic model to mimic a variety of ONA backbone chemistries, we do not precisely match the experimentally observed size and mass of DNA base pairs. We also adopt a generic CG model for polymer and ions such that the PL and IN beads do not represent a specific chemistry. To model the WC H-bonds, the complementary HB sites interact with the Lennard-Jones (LJ) potential

conjugation of ONA to polymers or small molecule linkers impacts ONA hybridization/melting. There are many studies on how DNA hybridization/melting are impacted by conjugation to other large/small molecules or grafting to substrates. Many experimental and theoretical studies show that conjugation of nucleic acids to nanoparticles52−54 and/or crosslinker molecules33,52,55−58 increases DNA stability and sharpens the DNA melting curve in the resulting self-assembled nanostructures. For the cross-linker molecules, the extent of this impact depends on the flexibility and architecture of crosslinkers. Improved stability and specificity of DNA duplex in these assembled structures was attributed to the high local density of DNA strands resulting in DNA cooperative binding. Multiple experimental59−62 and simulation63−70 studies show that conjugation with polymers, without forming a selfassembled structure, also affects the secondary structure of biomolecules, such as isolated peptides, and that such conjugation can favor or disfavor peptides’ secondary structure depending on the sequence and chemistry of peptides. To the best of our knowledge, there are no such studies aimed at elucidating effects of polymer conjugation on ONA hybridization for the different new ONA chemistries beyond the traditional DNA. In this work, we perform coarse-grained (CG) and all-atom (AA) simulations to understand the hybridization/melting thermodynamics of polymer−ONA conjugates in solution for varying ONA backbone flexibility and charge. Our CG approach allows us to explore experimentally relevant time and length scales with high computational efficiency. We also employ advanced sampling techniques, including parallel tempering71 and metadynamics simulations,72 to ensure efficient sampling of the configurational space. We also perform biased atomistic simulations to further appreciate the strengths and limitations of our CG model. In all approaches, we find that the 8−100-mer linear polymer conjugation further destabilizes the 8-mer ONA duplex in the case of less stable duplexes, while in stable duplexes the same polymer conjugation has minimal impact on ONA duplex stability. The paper is organized as follows. In section II, we describe the details of the CG model and the AA force field used in the CG and AA simulations, respectively. In section III, we present the details of the CG and AA simulations. In section IV, we explain our methods in analyzing the simulation data, and in section V, we summarize the parameters varied in this work. In section VI, we present the results of the CG and AA simulations, and then bridge the outcomes in the Discussion and Conclusion section at the end.

⎡⎛ HB ⎞12 ⎛ HB ⎞6 ⎤ σ σ UijHB ‐ complementary(r ) = 4εijHB⎢⎜ ⎟⎥ ⎟ −⎜ ⎢⎣⎝ r ⎠ ⎝ r ⎠ ⎥⎦

II. MODEL II.A. CG Model of ONA and Polymer. We recently developed a CG model to capture the hybridization/melting B

(1)


Article


Figure 2. (a) Chemical structure of PEG8-conjugated single-stranded DNA with sequence GCGC as used in AA-PT-MetaD-PT simulations. (b) The simulation snapshot of the resulting DNA duplex formed by two complementary PEG8-conjugated double-stranded DNA.

where εHB is the van der Waals energy parameter for the complementary HB sites. The subscript ij is dropped for σHB because all HB sites have the same size. The potential and force are both truncated at 1.9σHB. A GROMACS-style switch function74 is used to smoothly ramp potential and force to zero at 2σHB. We note that, despite the isotropic LJ potential between complementary HB sites, the directionality and specificity of WC H-bonds are captured by carefully designing relative size and position of BB and HB beads, and also by defining a repulsive shell that is only effective between HB sites of the same type.31 The interactions between all other pairs of beads (i.e., BB−BB, BB−HB, PL−PL, PL−BB, PL−HB, all interactions involving ions, and interactions between noncomplementary HB sites) are purely repulsive via Weeks− Chandler−Andersen (WCA) potential ⎧ ⎡⎛ σ ⎞12 ⎛ σ ⎞6 ⎤ ij ij 1/6 ⎪ ⎪ 4ε⎢⎝⎜ ⎠⎟ − ⎝⎜ ⎠⎟ ⎥ + ε r < σij(2) r r ⎦ Uij(r ) = ⎨ ⎣ ⎪ ⎪0 r ≥ σij(2)1/6 ⎩

model for water,78 and the Joung/Cheatham parametrization for ions.79 For the DNA moiety, we adopt the partial charges directly from Amber f f14sb. For the PEG moiety, we follow the proposed chemistry of Jeong et al.34 to attach an 8-mer of PEG (denoted as PEG8) to the 5′-terminal cytosine of DNA via an acid-cleavable phosphoramidate linkage (Figure 2). We adopt the general Amber force field80 (gaf f) for bonded and nonbonded interactions of the PEG8 and the linkage group, and then use the Antechamber package81 and Gaussian 09 software82 to calculate the partial charges (Supporting Information, section 1).

III. SIMULATION DETAILS In this work, we perform three types of simulations. In the first type, we perform unbiased coarse-grained simulations in a canonical ensemble (CG-NVT) using Langevin dynamics with the CG model described in section II.A. The CG-NVT simulation details resemble our recent publication,31 so we only describe the key features in section III.A. The advantages of CG-NVT simulations lie in the fact that they are relatively fast with reasonable computational expense and a straightforward simulation protocol, and allow us to conduct simulations of multiple ONA duplexes at the concentration of interest. This makes CG-NVT simulations a reasonable choice to explore our expanded parameter space. To ensure that the unbiased CG-NVT Langevin dynamics sufficiently samples the free energy landscape of the hybridization/melting process, we conduct a second type of simulation where we use CG Langevin dynamics with welltempered metadynamics83 (CG-NVT-MetaD) to ensure enhanced sampling of the configurational space along selected collective variables (CVs). In well-tempered metadynamics, Gaussians with time-dependent height are deposited to the configurational space to deter excessive sampling of the visited states (e.g., a local free energy minimum state). While details of the metadynamics approach are explained elsewhere,72 in section III.B, we only elaborate on our choice of CVs and other simulation related details. Although CG-NVT and CG-NVT-MetaD simulations take different approaches in sampling the configurational phase space, they are both bound to the inherent limitations of our CG model for ONAs. While we have verified our CG model against DNA melting curves in our recent publication,31 we still need to examine if and how excluding detailed polymer− ONA−water atomistic interactions impacts the results obtained using our CG model. Therefore, in the third type of simulations, we perform atomistic parallel tempering metadynamics simulations in the well-tempered ensemble21,84−86 (AAMetaD) to assess the effects of PEG conjugation and the

(2)

where the subscript ij is dropped for ε because all beads but the complementary HB sites have the same energy parameter (1.0ε). Arithmetic mean is used for unlike LJ size parameters. LJ interactions are also effective between all 1−3 and 1−4 bonded pairs to avoid overlap of beads in flexible chains. The purely repulsive interactions between all but complementary HB sites represent a system where the implicit solvent is equally good for both polymer and ONAs, and the polymer and ONA do not have any specific attractive interactions. Electrostatic interactions between charged beads (BB, IN) are calculated by the Coulomb potential qiqj Uijcoul(r ) = 4πε0rεr(T ) (3) where qi is the charged valency on bead i, ε0 is the permittivity of a vacuum, T is the temperature, and εr is the relative permittivity of implicit solvent, which at dilute salt concentrations is estimated via75 εr(T [K ]) = 249.4 − 0.788T + 0.00072T 2

(4)

where the dielectric constant is assumed to be distanceindependent and the temperature is in units of Kelvin. In our CG model, 1−3 and 1−4 electrostatic interactions between BB beads are set to zero. II.B. Atomistic Force Field of DNA and PEG. For atomistic simulations of DNA, we use the Amber f f14sb combination of force fields,76 which includes parmbsc0 modification for the nucleic acid backbone,77 the TIP3P C


Article


second CV is the number of HB−HB bonds between the two hybridizing strands, which is calculated using the following switch function

presence of explicit water on DNA hybridization/melting. While metadynamics enhances sampling of the configurational space only along selected CVs, parallel tempering enhances sampling along all degrees of freedom by exchanging configurations from higher temperatures. Performing AAMetaD simulations for experimentally relevant time and length scales, and duplex concentrations, is not feasible. Therefore, AA-MetaD simulations of this work, with details in section III.C, are presented here only as an auxiliary tool to understand the strengths and limitations of our CG model for future simulation studies of polymer−ONA based systems. III.A. CG-NVT Simulations. We perform CG molecular dynamics simulations in an NVT ensemble using the LAMMPS87 package. The initial configuration for each simulation is prepared by randomly placing fully hybridized duplexes in a cubic simulation boxwith periodic boundary conditionswhile avoiding overlaps. The number of duplexes and the simulation box size are chosen on the basis of the duplex concentration of interest. For simulations with 0.006 mM duplex concentration, 20 duplexes are randomly placed in a box of length 300σ. For simulations with 1.5 mM duplex concentration, 65 duplexes are randomly placed in a box of length 70σ. When BB beads are neutral, simulations are performed in the absence of salt. When BB beads are (negatively) charged, in addition to neutralizing positively charged counterions, 1 mM of monovalent salt is also added. These ions are randomly placed in the simulation box while avoiding particle overlaps. The initial configurations are subjected to Langevin dynamicş and the equations of motion are integrated using a two-level rRESPA multiple time step algorithm,88 with a time step of 0.001 and 0.0005 in reduced units (corresponding to 6 and 3 fs, respectively) for nonbonded and bonded interactions, respectively. The temperature is controlled by setting the friction coefficient to 10 (in reduced units of time). Electrostatic interactions are treated with the particle− particle−particle−mesh (PPPM) method,89 with a force tolerance of 10−2, interpolation order of 2, and real space cutoff distance of 20σ. The equilibrium runs vary between 1.0 × 108 and 1.4 × 108 time steps, corresponding to 600−840 ns. The production runs vary between 1.0 × 107 and 2.0 × 107 time steps, depending on the system under study. Longer equilibrium and production times are needed for simulations around the melting temperature of each system. In the production runs, configurations are stored every 2000 time steps for data analyses. Three independent simulation trials, with different initial velocities and random number seeds, are performed to estimate statistical uncertainties. III.B. CG-NVT-MetaD Simulations. In the CG-NVTMetaD simulations, we simulate only one duplex in a cubic simulation box of length 50σ with periodic boundary conditions. The simulation details of CG-NVT-MetaD are similar to CG-NVT simulations, but the initial configurations are relaxed only for 105 time steps, and then Gaussians are deposited along two CVs using the Plumed 2.2 plugin.90 The first CV, denoted by δr, measures the separation between the hybridizing ONA single strands, and is the distance between the centers of mass (COM) of each strand. To calculate the COM, the terminal bases are excluded due to their excessive fluctuations. A hypothetical wall is placed at δr = 15σ to avoid spending simulation time sampling fully melted states. The wall is a harmonic potential with a force constant of 1500ε/σ2, and is effective whenever the first CV is greater than 15σ. The

N HB =

1 − (r /r0)n 1 − (r /r0)m

(5)

where NHB is the HB−HB coordination number, r is the distance between any two HB sites in a duplex, n = 100, m = 120, and r0 = 0.5σ. This differentiable form of this switch function gradually allows conformations to move toward unexplored values of this reaction coordinate. The parameters of eq 5 determine the shape of the switching function, and are chosen so that NHB = 8.00 and 0.00 for fully hybridized and fully melted 8-mer duplexes, respectively. In the framework of our CG model, the two CVs above are sufficient to distinguish every possible ONA conformation during the hybridization/melting process. The bias factor in metadynamics simulations is set to 5, and the Gaussians are deposited every 1000 time steps. The initial height and width of Gaussians for the first and second CVs are 10ε and 0.1 (in units of CV), respectively. The production run is continued for about 6 × 108 time steps until the free energy landscape converges. The convergence is assessed by monitoring the free energy of hybridization over simulation time. Three independent trials are performed, with different initial velocities and random number seeds for the Langevin dynamics, to estimate statistical uncertainties. The free energy landscape obtained with these single duplex CG-NVT-MetaD simulations is corrected to estimate the landscape for the duplex concentration of interest. This is described in the section IV. III.C. AA-MetaD Simulations. We use Gromacs 5.191 and the Plumed 2.2 plugin90 to run AA-MetaD simulations. Using the Avogadro molecular editor,92 two separate sets of initial configurations are prepared: B-DNA duplex with sequence GCGC and (double-stranded) B-DNA-PEG8 conjugates with sequence GCGC-PEG8. The following steps are performed on each system. First, the initial configuration is solvated in a cubic simulation box of length 6.42 nm, and counterions are added, using the internal functionalities of the Gromacs package. Second, the initial configurations are relaxed using multiple short simulations in NPT and NVT ensembles (Supporting Information, section 2). Third, six replicas of each system are prepared at temperatures of 300, 316, 334, 354, 376, and 400 K. Replicas are further relaxed at the corresponding temperatures by running a 0.2 ns parallel tempering simulation in the NVT ensemble. At this stage, due to the large temperature intervals, no exchange of configurations occurs between the neighboring replicas. Fourth, a 1 ns well-tempered metadynamics (WTM) simulation is performed using the following collective variables: the total number of WC H-bonds between the two hybridizing strands (NWC H‑bonds) and the total number of stacking interactions in the duplex (NStacking). The NWC H‑bonds is measured using the switching function of eq 5, with n = 8, m = 12, and r0 = 2.5 Å, where r is the distance between a WCrelevant hydrogen atom in one strand and the corresponding heavy atom in the opposite strand.21 While calculating NWC H‑bonds, all possible base pair hybridizations are considered; i.e., a cytosine base in the first strand can hybridize with either guanine base in the second strand. For computational efficiency, the complete melting of the duplex is prevented by placing a hypothetical wall (harmonic potential, with a force constant of 1000 kj/mol/nm2) at NWC H‑bonds = 1. To quantify D


Article


A(δr , N HB) = Aedge(δr , N HB) + 2kBT ln(δredge) − 2kBT ln(δr )

the second CV, we use the functional form of eq 5, with n = 6, m = 12, and r0 = 6 Å, where r is the distance between the glycosidic nitrogen atom in neighboring bases on one strand. In a 4-mer duplex, each strand has a maximum of 6 stacking interactions (by double-counting the stacking of the central CG), which leads to a maximum of 12 stacking interactions for the entire DNA duplex. This short WTM simulation encourages the replicas to diffuse along the aforementioned CVs, and generate unique initial states for each replica in the CV phase space.93 The accumulated bias along the CVs at this step is not used in the subsequent steps. Fifth, a 3 ns parallel tempering metadynamics simulation is performed using the total potential energy as CV. This step perturbs the potential energy of the replicas, and provides enough overlaps in the potential energies of the neighboring temperatures so that the exchange acceptance rate between neighboring replicas boosts to 30−40%. This concludes the equilibration stage for the AAMetaD simulations. Lastly, using the potential energy of the last step as a static bias and the CVs defined for the fourth step, AAMetaD simulations are performed on each system for 500 ns. The metadynamics parameters are given in the Supporting Information, section 3. We note that, even though the parallel tempering is applied in all the steps above, the exchanges between temperatures are only successful after perturbing the potential energy. For each system, the convergence of the free energy is assessed by examining the relative depths of the free energy minima, similar to a recent publication.21 Two independent trials are performed for each system to evaluate statistical uncertainties.

(6)

where A is the Helmholtz free energy, T is the temperature, kB is the Boltzmann number, and Aedge is the free energy at δredge. Having obtained the free energy landscape at the concentration of interest (same as the CG-NVT simulations), we extract the duplex melting curve. We do this by first placing a saddle line at NHB = 4 (which is half of the maximum number of H-bonds in an 8-mer duplex) and calculating the probabilities of hybridized and melted states using the Boltzmann factor Pj =

∑ exp(−Aj /kBT );

j ∈ hybridized, melted (7)

n

where P is the probability and the summation runs over either hybridized or melted states. Once the probabilities of hybridized and melted states are known, f50% is calculated by using the method of Ouldridge et al.95 f50% = 1 +

1 − 2Φ

⎛ 1 ⎞⎟2 ⎜1 + −1 ⎝ 2Φ ⎠

(8)

where Φ is the ratio of the probabilities of hybridized and melted states (Phybridized/Pmelted). From the probabilities calculated in eq 7, the free energy change upon hybridization is ΔAhybridization = kBT ln(P hybridized /P melted)

(9)

which is simply the sum over the free energy grid for the hybridized and melted regions of the free energy landscape. In the AA-MetaD simulations, we only consider the configurations and the hill file of the lowest temperature (300 K). The configurations at higher temperatures are only used to enhance sampling of the configurational space in the basin with the lowest temperature, and thus are not used in the analyses. The two-dimensional free energy landscape is generated by building histograms along the two CVs (NWC H‑bonds, Nstacking) using a histogram bin size of 0.2 for both CVs. To assess the convergence of PT-MetaD-WT simulations, one-dimensional free energies are generated by integrating along Nstacking, and the free energy differences between relevant local minima are monitored over simulation time. We assume the AA-metaD simulations are converged if the differences between free energy local minima reach a plateau, and the standard deviation over the two independent trials is not greater than 1 kcal/mol. To generate the two-dimensional free energy for each nucleobase, the histogram reweighting approach of Bonomi et al.96 is used to postprocess the hill file of the lowest temperature. The definitions of CVs for each nucleobase are the same as those explained in section III.C but are defined for a nucleobase rather than the entire duplex. The histograms are calculated by setting the bin sizes of 0.04 and 0.02 and band widths of 0.08 and 0.04 for NWC H‑bonds and Nstacking, respectively. The resulting free energies are summed separately for the terminal and central bases. To estimate the free energy of WC H-bonding (ΔAWC H‑bonds) and stacking (ΔAstacking) for the central and terminal bases, eq 7 is used to calculate the relevant probabilities, and then, eq 9 is used to calculate the free energy difference between bound and unbound as well as stacked and unstacked states. The saddle line in each case is defined at the center of the CV range. For example, each G···C base pair forms three WC H-bonds, and the saddle line is placed at NWC H‑bonds = 1.5 for calculating ΔAWC H−bonds. In the

IV. ANALYSES In the CG-NVT simulations, the analyses for calculating ONA melting curves and H-bond residence time distributions are the same as those described in our recent publication.31 Melting curves are obtained by obtaining an ensemble average fraction of duplexes in the simulation box that are hybridized 50% or more (f50%). The melting temperature (Tm) of an ONA is the temperature at which f50% = 0.5. The H-bond residence time plot serves as a proxy for the extent of the base-pair breathing, and is estimated by calculating the ensemble average fraction of configurations where the base pairs in the hybridized duplexes maintain their H-bonds (τhybridized). In the CG-NVT-MetaD simulations, we use the sum_hills tool of the Plumed plugin to calculate the two-dimensional histogram of the free energy along the two CVs described in section III.B. To calculate the histogram, the bin sizes are set to 0.1σ and 0.3 for δr and NHB, respectively. As described in section III.B, we run the CG-NVT-MetaD simulations with only one duplex, and place a hypothetical wall at δr = 15σ to increase computational efficiency. Because this setup corresponds to a higher duplex concentration than those studied by CG-NVT simulations, we need to extend the free energy landscape to distances beyond δr = 15σ to retrieve the same concentrations as CG-NVT simulations. To extend the free energy landscape beyond δr = 15σ, we realize that the other CV for CG-NVT-MetaD simulations, NHB, is zero for δr > 15σ, meaning that, beyond the hypothetical wall, the free energy is only a function of δr, and is dominated only by the entropy of the strands. Therefore, following de Pablo and co-workers94 for extending the free energy beyond δr = δredge = 15σ, the free energy expression becomes E


Article


method, and bridge the results of different simulation types together. VI.A. Effects of Polymer Conjugation Studied by CGNVT Simulations. We first assess the effects of polymer conjugation on the melting curve of ONAs with a charged sflx ONA backbone. For 8-mer charged sflx ONA duplexes with GCGCGCGC sequence, conjugation with (8−40-mer) linear polymers and 15-mer 4-arm star-like polymers with varying flexibilities studied here exhibits no significant impact on the melting temperature and the melting transition width (Figure 3a). However, for 8-mer charged sflx ONA duplexes with

analyses for the central and terminal bases, the configurations at which the duplex is in the searching mode are excluded, where the searching mode is defined as a duplex state at which the two hybridizing strands slide along the duplex. The distributions of end-to-end distance (Ree) and the radius of gyration (Rg) for the PEG chain are also calculated to understand the effects of ONA conjugation on the conformation of the polymer moiety. Ree is defined as the distance between C1 and O1 atoms in the PEG molecular structure (refer to Figure S1 in the Supporting Information), and Rg is calculated on the basis of all atoms in the PEG moiety. The bin sizes for calculating the Ree and Rg histograms are 0.06 and 0.02 nm, respectively.

V. PARAMETERS VARIED IN THIS PAPER In the CG-NVT simulations, we vary the conjugated polymer length, architecture, and flexibility and ONA sequence, ONA backbone charge, and ONA flexibility. To limit the parameter space, we limit the simulations to two ONA sequences, GCGCGCGC and ATATATAT, which represent strands with 100 and 0% G···C contents, respectively. We also focus only on two ONA backbone types: (1) charged and semiflexible backbones (kangle = 30 ε/rad2) denoted as sflx; (2) neutral and flexible backbones (kangle = 10 ε/rad2) denoted as f lx. The former mimics a DNA-like backbone, and the latter mimics a PNA73 (or CNA45)-like backbone. We consider two duplex concentrations (c), 0.0058 and 1.5 mM, which cover a wide range of practical applications of ONAs in solutions. For the polymer moiety, we change the polymer length from 8 to 100 beads, polymer flexibility from fully flexible (kangle = 10 ε/rad2), denoted as f f lx, to semiflexible (kangle = 30 ε/rad2), denoted as sflx, and also change the polymer architecture from linear to star with 4-arms, each arm being a 15-mer f flx polymer. The polymers are attached or conjugated to ONA, as shown in Figure 1. In CG-NVT-MetaD simulations, we only focus on a few systems to verify the true sampling of the configurational space by CG-NVT simulation. Performing CG-NVT-MetaD simulations for all of the systems that we studied by CG-NVT simulations is neither computationally efficient nor necessary. The studied systems by CG-NVT-MetaD simulations are (1) neutral f lx ONA with sequence CGCGCGCG at c = 0.0058 mM, (2) the same as in 1 but with sflx ONA backbone, (3) the same as in 2 but with the nucleobase distance being 1.2 times larger than that of DNA, (4) neutral flx ONA with sequence CGCGCGCG at c = 0.0058 mM conjugated to f f lx linear polymers with length 8, and (5) the same as in 4 but with the polymer length being 40. In AA-MetaD simulations, we have two systems: (1) 4-mer DNA duplex with sequence CGCG; (2) the same as in 1 but conjugation with an 8-mer PEG according to chemistry shown in Figure 2. We also perform simple AA-NVT simulations on PEG8 without using metadynamics and parallel tempering biasing approaches to assess the effects of conjugation on polymer conformation.

Figure 3. Effects of polymer length, flexibility, and architecture on the melting curves of charged sf lx ONAs with (a) GCGCGCGC and (b) ATATATAT sequences at c = 0.0058 mM obtained from CG-NVT simulations. Black squares (with fitted lines) represent the melting curves for duplexes with no conjugated polymer. The melting curves of polymer-conjugated duplexes are shown with blue circles (8-mer sflx linear polymer), green upward triangles (8-mer f f lx linear polymer), red downward triangles (40-mer sflx linear polymer), brown diamonds (40-mer f f lx linear polymer), and light blue stars (15-mer 4-arm f f lx star polymer). The abbreviations sflx and f f lx stand for semiflexible and fully flexible backbones, respectively. The solid and dashed vertical lines in part b mark the Tm for duplexes with no polymer and 15-mer 4-arm f f lx star polymer, respectively.

ATATATAT sequence, conjugation with 8-mer and 40-mer f f lx linear polymers, and 15-mer ff lx star-like polymer, slightly destabilizes the duplexes marked by a drop in Tm by ∼2.5 ± 1 °C (Figure 3b). We emphasize that the difference in Tm (shown with dashed lines in Figure 3b) is obtained from the melting temperatures calculated using the entire melting curves rather than just one temperature. Our reasoning for the results observed in Figure 3 is as follows: The conjugation of ONAs with purely repulsive polymers alters the change in conformational entropy upon hybridization. Because the gain in enthalpy upon hybridization is lower for the ATATATAT sequence than the GCGCGCGC sequence, the contribution of entropy to the free energy change upon hybridization (ΔAhybridization) becomes more significant for the ATATATAT sequence than the GCGCGCGC sequence. Therefore, the additional contribution to free energy coming from the loss in conformational entropy of polymer upon duplex hybridization appears to have a more pronounced effect on ΔAhybridization of the ATATATAT sequence. Furthermore, at higher duplex concentration due to the relatively smaller available free volume for duplexes in the system, the loss in conformational entropy is smaller, so the conjugation with polymer imposes an insignificant impact on the melting curves for both ATATATAT and GCGCGCGC sequences (Supporting Information, Figure S2). Lastly, the 8− 40-mer linear and 15-mer 4-arm star athermal polymers all have sizes (e.g., ideal chain radius of gyration or end−end distances) that are commensurate with the 8-mer ONA duplexes.

VI. RESULTS In this section, we present the results from the CG-NVT, CGNVT-MetaD, and AA-MetaD simulations for the effects of polymer conjugation on ONA hybridization/melting. We also discuss the limitations and strengths of each simulation F


Article

The Journal of Physical Chemistry B Therefore, these systems do not show variations in the effects described above for varying polymer lengths or architectures considered here. Our CG-NVT simulations for ONAs with neutral f lx backbones show that conjugation with the polymers considered in this study has no significant impact on the melting curves and the melting transition widths of duplexes with GCGCGCGC (Figure 4a) and ATATATAT (Figure 4b)

Figure 5. Effects of polymer conjugation on the fraction of time that base pairs remain hybridized (τhybridized) in hybridized, neutral f lx ONA duplexes with GCGCGCGC sequence at c = 0.0058 mM. Black squares (with fitted lines) represent the data for duplexes with no conjugated polymer. The data for polymer-conjugated duplexes are shown with green upward triangles (8-mer f f lx linear polymer), brown diamonds (40-mer f flx linear polymer), and dark green pentagons (100-mer f f lx linear polymer).

Figure 4. Same as Figure 3 but for ONAs with neutral flx backbones. The dark green pentagon symbols in part a show simulation data for ONA conjugated to 100-mer f flx linear polymer.

sequences at c = 0.0058 mM. The same statement also holds for c = 1.5 mM (Supporting Information, Figure S3). For duplexes with neutral backbones, the absence of electrostatic repulsion between hybridizing strands leads to significantly more favorable enthalpy change upon hybridization compared to the same sequences with the charged backbones. Therefore, we conjecture that the energetic contribution to ΔAhybridization dominates over the entropic contribution, so that the additional loss in conformational entropy brought about by polymer conjugation has no effect on the hybridization thermodynamics of either sequences that we studied in this work. Despite the small to insignificant impact of polymer conjugation on the ONA melting temperature, our analyses on GCGCGCGC duplexes with neutral f lx (Figure 5) and charged sf lx ONA backbones (Supporting Information, Figure S4) show that conjugation with polymer decreases the extent of the base-pair breathing by stabilizing the hydrogen bonds between complementary bases in the hybridized duplexes. This is more pronounced for neutral flx ONA backbones, especially in temperatures above the melting temperatures of the studied duplexes. As stated earlier, there is an additional entropic contribution, the loss in conformational entropy of the polymer upon hybridization, that has to be compensated with the energetic gain upon hybridization. Any base breathing or “bubbles” in the duplex reduces that energetic gain. Therefore, the base breathing within the duplex is reduced upon polymer conjugation. However, this stabilizing effect is insufficient to significantly alter the Tm values of ONAs upon conjugation with polymers. The same conclusion also holds for duplexes with ATATATAT sequence (Supporting Information, Figure S5). VI.B. Effects of Polymer Conjugation Studied by CGNVT-MetaD Simulations. To ensure that the results presented in section VI.A are not impacted by insufficient sampling of an unbiased molecular dynamics simulation, we investigate the effect of polymer conjugation on ONA duplex stability for a select few systems with NVT-MetaD simulations using the same CG model, as described in the Simulation

Details section. First, we assess the convergence of the CGNVT-MetaD simulations (Supporting Information, Figure S6) as we reproduce some of the conclusions that we obtained in our recent publication31 for the effects of ONA backbone flexibility and nucleobase distance on the ONA hybridization/ melting curves in the absence of polymer conjugation. We demonstrate that, in the absence of polymer conjugation, the CG-NVT-MetaD simulations reproduce our published CGNVT simulations results which show that for 8-mer ONA with GCGCGCGC sequence, changing backbone flexibility from f lx to sflx and nucleobase distance from 0.84σ to 1.0σ increases Tm by ∼40 and ∼5 °C, respectively (Supporting Information, Figure S7). These comparisons confirm that at least for the ONA systems without polymer conjugation the CG-NVT simulations and CG-NVT-MetaD simulations are consistent with each other. We then perform CG-NVT-MetaD simulations for selected systems presented in section VI.A to capture the effects of polymer conjugation on hybridization/ melting thermodynamics. We specifically study how the melting curves of neutral f lx ONAs with GCGCGCGC sequence change upon conjugation with 8-mer, and 40-mer, f f lx linear polymers (Figure 6). The simulations show that the free energy landscape has two local minima at (dr ≈ 1σ, NHB > ∼6) and (dr > ∼4σ, NHB ≈ 0), which correspond to the fully hybridized and fully melted duplexes, respectively. The relative depth of these minima changes with temperature, and determines which duplex state is more favorable. Our analysis on these free energy landscapes shows that conjugation with 8-mer and 40mer f f lx linear polymers changes Tm by approximately −1 and +1 °C, respectively, compared to the duplex with no polymer attached. However, these small changes are within the statistical uncertainty of our estimation for Tm based on three independent simulation trials. Therefore, we conclude that CG-NVT-MetaD simulations are consistent with CG-NVT simulations in showing that conjugation with polymers has no significant impact on the hybridization/melting of duplexes with GCGCGCGC sequence. G


Article


Figure 6. Effects of polymer conjugation on the free energy landscape (a, b) and the melting curves of ONA duplexes with a neutral f lx backbone and GCGCGCGC sequence obtained from CG-NVT-MetaD simulations. (a) The free energy landscape at T = 54 °C for the ONA duplex with no polymer attached. The inset shows representative snapshots of the fully hybridized and the fully melted duplexes, and the approximate regions of the free energy landscape that corresponds to these snapshots. The dashed line illustrates the saddle line used in eqs 7 and 8 to distinguish between hybridized and melted states. (b) The free energy landscape at T = 54 °C for the ONA duplex conjugated with 8-mer f f lx linear polymer. (c) The melting curves obtained from GC-NVT-MetaD simulations. Black squares (with fitted lines) represent the melting curve for the ONA duplex with no polymer attached. The melting curves of the polymer-conjugated duplexes are shown with green upward triangles (8-mer f f lx linear polymer) and brown diamonds (40-mer f f lx linear polymer). The gold star marks the data points for which the free energy landscapes are shown.

Figure 7. Effects of polymer conjugation on the free energy landscape of the central bases in the GCGC (a) and GCGC-PEG8 (b) systems. The free energy is averaged over the central bases that are marked with green dashed ovals in the representative simulation snapshot for each system. In the snapshots, cytosine, guanine, and PEG residues are colored with red, blue, and gray, respectively, for further clarity. The local minima (1), (2), and (3) in the free energy landscape of part a correspond to different possible DNA conformations, few of which are depicted with the cartoon representation of DNA duplex on the right column.

visual inspection of the free energy landscape shows the configurational states at which the DNA duplex is melted are visited with higher probability in the GCGC−PEG8 than the GCGC system. To quantify this observation, we extract onedimensional free energies as a function of NWC H‑bonds and Nstacking from the two-dimensional free energy landscape (Supporting Information, Figure S10). For the GCGC duplex, there are three distinct minima at NWC H‑bonds ≈ 6, 9, and 11 and the stacking interactions are more favorable in the Nstacking ≈ 10 region. For the GCGC−PEG8 duplex, the three minima are at NWC H‑bonds ≈ 1, 3, and 6 and stacking interactions are more favorable in the Nstacking < 6 region. The shift to lower values of NWC H‑bonds and Nstacking upon PEG conjugation suggests that conjugation with PEG destabilizes the DNA duplex. Due to the expected excessive fraying of the terminal bases, especially in a short 4-mer DNA duplex, we focus on the effects

While both CG-NVT-MetaD simulations and CG-NVT simulation methods show consistency in our conclusion, they are based on the same CG model. While this CG model correctly captured the melting curves seen for DNA duplexes without polymer conjugation in our recent publication,31 it is important to understand the limitations of the CG model in capturing the polymer conjugation effect. To achieve this, we conduct atomistic metadynamics simulations with DNA (as the atomistic force fields are well tested for this oligonucleic acid) with explicit water and choice of PEG for the conjugated polymer. VI.C. Effects of PEG Polymer Conjugation Studied by AA-MetaD Simulations. Upon convergence of the relative depths of the free energy minima over time (Supporting Information, Figure S8), we compare the free energy landscape with the collective variables NWC H‑bonds and Nstacking in duplexes with and without PEG8 (Supporting Information, Figure S9). A H


Article


temperatures of DNA with GCGC and ATATATAT sequences are around −100 and −20 °C, respectively97) conjugation with polymers of sizes commensurate with the duplex further destabilizes the duplex. On the other hand, this seemingly consistent conclusion seems to stem from different molecularlevel reasoning in each of these methods. In our coarse-grained approach, we model generic polymers with purely repulsive interactions, so the impact on the stability of ONAs with ATATATAT sequence is attributed to the additional loss in conformational entropy of the polymer upon hybridization that becomes significant only compared to the weak energetic gain upon hybridization of an ATATATAT duplex. In the CG model, the solvent is not modeled explicitly and any entropic/ energetic contributions from the direct solvent−ONA or solvent−polymer interactions are not captured. In contrast, in the atomistic simulations, we study a specific chemistry of polymer, PEG, where we find that the PEG moiety collapses on the DNA duplex, shielding the relative hydrophobic nucleobases from water, and increasing the effective electrostatic repulsion between the hybridizing DNA strands, which consequently destabilizes the DNA duplex with GCGC sequence. Furthermore, unlike the coarse-grained simulations, in the atomistic simulations, we are limited to a short PEG chain and a short DNA duplex (4 base pairs), while the coarsegrained simulations explore a longer DNA duplex (8 base pairs) and longer (generic) polymer chains (8-mer to 40-mer) and star and linear architectures. The coarse-grained generic polymer−implicit solvent pair is not a good mimic of PEG− water interactions.70,98 Keeping these differences between the two approaches in mind, we conclude that in experimental systems both the conformational entropy loss of the polymer as well as the direct interactions of the ONA−polymer conjugate with the solvent will dictate the thermal stability of the polymer conjugated ONA duplex. Coarse-grained and atomistic simulations are both bound to several limitations when used for studying the effects of polymer conjugation on the ONA hybridization, and therefore, care must be taken before generalizing their outcomes. The atomistic simulations in explicit water provide a detailed understanding of DNA/water, DNA/PEG, and PEG/water interactions but suffer from (1) low computational efficiency, which in our case forces us to simulate one 4-mer DNA duplex; (2) uncertainty in the viability of the general Amber (or other atomistic) force field for DNA−PEG conjugates, which is particularly important in a 4-mer DNA because half of the nucleobases are directly attached to the PEG moiety; (3) uncertainty in relating the extent of base-pair breathing (as studied in this work via monitoring nucleobase stacking and Hbonding interactions) to the melting curve and the melting temperature of the duplex. Due to these limitations, it is not straightforward to extrapolate the conclusions of atomistic simulations with 4-mer DNA duplex to experimentally relevant systems, where longer duplexes are used. In contrast, the coarse-grained simulations in implicit solvent are computationally efficient, enabling studies of longer duplexes and experimentally relevant concentrations, but are limited to inherent drawbacks of our coarse-grained model for polymerconjugated ONAs, which are (1) coarse-graining each nucleobase into one backbone (BB) and one H-bond (HB) bead, and thus missing detailed atomistic interactions; (2) not representing the true (and fairly complex) PEG behavior70,98 by simple coarse-graining of the polymer moiety with the purely repulsive interactions in implicit solvent; (3) absence of explicit

of PEG conjugation purely on the central bases, rather than the entire duplex. Figure 7 shows the free energy landscape for only the central bases in DNA (a) and PEG8−DNA (b) duplexes. For both duplexes, there are three local minima at (NWC H‑bonds ≈ 3, Nstacking ≈ 2), (NWC H‑bonds ≈ 0, Nstacking ≈ 1), and (NWC H‑bonds ≈ 0, Nstacking ≈ 0). To quantify the effects of PEG conjugation, we calculate ΔAWC H‑bonds and ΔAstacking, according to the procedure explained in section IV. Our analyses show that ΔAWC H‑bonds for the central bases in DNA and PEG8− DNA duplexes is −0.36 ± 0.10 and 0.24 ± 0.12 kcal/mol, respectively. This means conjugation with PEG disfavors formation of WC H-bonds between complementary central bases. Similarly, ΔAstacking for the central bases increases from −0.52 ± 0.11 kcal/mol in DNA to 0.15 ± 0.06 kcal/mol in DNA−PEG duplexes. This also shows that stacking interactions become less favorable in PEG8−DNA compared to DNA. Because WC H-bonding and stacking are the two main driving forces for hybridization, we conclude that conjugation with PEG8 destabilizes the 4-mer DNA duplex that is expected to be unstable at 300 K. To understand what drives the conjugation with PEG8 to destabilize the 4-mer DNA duplex at 300 K, we analyze the PEG conformation via the Ree distribution; it shows a shift toward smaller values of Ree in DNA−PEG duplex compared to the free PEG in water (Supporting Information, Figure S11). The radial distribution function of terminal guanine and water shows that the total number of guanine−water pairs decreases in DNA−PEG duplex compared to the DNA duplex (Supporting Information, Figure S12). Both the PEG conformation and guanine−water pair correlation results suggest that the PEG chains collapse on the DNA, which is likely to shield the relatively hydrophobic DNA nucleobases from water. In other words, the PEG chains are likely solubilizing the DNA. Fewer number of water molecules in the vicinity of the DNA bases leads to a lower effective dielectric constant for solution, which in turn increases the effective electrostatic repulsion between the nucleobases in the hybridizing DNA strands. Increasing electrostatic repulsion between the nucleobases should destabilize the base pair hybridization for the PEG−DNA duplex compared to the DNA duplex.

VII. DISCUSSION AND CONCLUSION In this work, using coarse-grained simulations, we found that conjugation with purely repulsive 8−100-mer polymers destabilizes 8-mer oligonucleic acid duplexes with weaker WC H-bonding interactions at low duplex concentrations, while it has an insignificant impact on duplexes with stronger WC Hbonding. This conclusion holds for a range of polymers with varying flexibility and lengths (8-mer to 40-mers) and architecture (linear and four-arm star) studied in this paper. Atomistic simulations of a 4-mer DNA duplex in explicit water with and without conjugation to 8-mer PEG showed that for a short DNA duplex at T = 300 K (where it is expected to be weakly held together at this temperature) the conjugation with PEG destabilizes DNA duplex. These results are valuable for scientists who synthesize and create supramolecular assemblies/gels of oligomers/polymers attached to oligonucleic acids for engineering artificial tissues, drug or DNA delivery vehicles, and devices for bionanotechnology applications. On the one hand, atomistic and coarse-grained simulations are consistent in that they show for duplexes with weak enthalpic driving force of hybridization (the predicted melting I


Article


(13) Suzuki, Y.; Endo, M.; Sugiyama, H. Mimicking MembraneRelated Biological Events by DNA Origami Nanotechnology. ACS Nano 2015, 9, 3418. (14) Kocabey, S.; Kempter, S.; List, J.; Xing, Y.; Bae, W.; Schiffels, D.; Shih, W. M.; Simmel, F. C.; Liedl, T. Membrane-Assisted Growth of DNA Origami Nanostructure Arrays. ACS Nano 2015, 9, 3530. (15) Breslauer, K. J.; Frank, R.; Blöcker, H.; Marky, L. A. Predicting DNA duplex stability from the base sequence. Proc. Natl. Acad. Sci. U. S. A. 1986, 83 (11), 3746−3750. (16) Sugimoto, N.; Nakano, S.-i.; Yoneyama, M.; Honda, K.-i. Improved thermodynamic parameters and helix initiation factor to predict stability of DNA duplexes. Nucleic Acids Res. 1996, 24 (22), 4501−4505. (17) SantaLucia, J. A unified view of polymer, dumbbell, and oligonucleotide DNA nearest-neighbor thermodynamics. Proc. Natl. Acad. Sci. U. S. A. 1998, 95 (4), 1460−1465. (18) Sugimoto, N.; Satoh, N.; Yasuda, K.; Nakano, S.-i. Stabilization factors affecting duplex formation of peptide nucleic acid with DNA. Biochemistry 2001, 40 (29), 8444−8451. (19) Owczarzy, R.; You, Y.; Groth, C. L.; Tataurov, A. V. Stability and Mismatch Discrimination of Locked Nucleic Acid−DNA Duplexes. Biochemistry 2011, 50 (43), 9352−9367. (20) Hudson, G. A.; Bloomingdale, R. J.; Znosko, B. M. Thermodynamic contribution and nearest-neighbor parameters of pseudouridine-adenosine base pairs in oligoribonucleotides. RNA 2013, 19 (11), 1474−1482. (21) Elder, R. M.; Pfaendtner, J.; Jayaraman, A. Effect of hydrophobic and hydrophilic surfaces on the stability of double-stranded DNA. Biomacromolecules 2015, 16 (6), 1862−1869. (22) Singh, Y.; Murat, P.; Defrancq, E. Recent developments in oligonucleotide conjugation. Chem. Soc. Rev. 2010, 39 (6), 2054−2070. (23) Pinheiro, V. B.; Holliger, P. The XNA world: progress towards replication and evolution of synthetic genetic polymers. Curr. Opin. Chem. Biol. 2012, 16 (3), 245−252. (24) Pinheiro, V. B.; Taylor, A. I.; Cozens, C.; Abramov, M.; Renders, M.; Zhang, S.; Chaput, J. C.; Wengel, J.; Peak-Chew, S.-Y.; McLaughlin, S. H. Synthetic genetic polymers capable of heredity and evolution. Science 2012, 336 (6079), 341−344. (25) Pinheiro, V. B.; Holliger, P. Towards XNA nanotechnology: new materials from synthetic genetic polymers. Trends Biotechnol. 2014, 32 (6), 321−328. (26) Taylor, A. I.; Arangundy-Franklin, S.; Holliger, P. Towards applications of synthetic genetic polymers in diagnosis and therapy. Curr. Opin. Chem. Biol. 2014, 22, 79−84. (27) Takezawa, Y.; Shionoya, M. Metal-mediated DNA base pairing: alternatives to hydrogen-bonded Watson−Crick base pairs. Acc. Chem. Res. 2012, 45 (12), 2066−2076. (28) Hirao, I.; Kimoto, M.; Yamashige, R. Natural versus artificial creation of base pairs in DNA: origin of nucleobases from the perspectives of unnatural base pair studies. Acc. Chem. Res. 2012, 45 (12), 2055−2065. (29) Harcourt, E. M.; Kool, E. T. Designer bases, base pairs, and genetic sets: biochemical and biological activity. Synthetic Biology; Royal Society of Chemistry: 2014; Vol. 1, pp 1−30. (30) Dhami, K.; Malyshev, D. A.; Ordoukhanian, P.; Kubelka, T.; Hocek, M.; Romesberg, F. E. Systematic exploration of a class of hydrophobic unnatural base pairs yields multiple new candidates for the expansion of the genetic alphabet. Nucleic Acids Res. 2014, 42 (16), 10235−10244. (31) Ghobadi, A. F.; Jayaraman, A. Effect of backbone chemistry on hybridization thermodynamics of oligonucleic acids: a coarse-grained molecular dynamics simulation study. Soft Matter 2016, 12, 2276− 2287. (32) Immoos, C. E.; Lee, S. J.; Grinstaff, M. W. DNA-PEG-DNA triblock macromolecules for reagentless DNA detection. J. Am. Chem. Soc. 2004, 126 (35), 10814−10815. (33) Gibbs, J. M.; Park, S.-J.; Anderson, D. R.; Watson, K. J.; Mirkin, C. A.; Nguyen, S. T. Polymer-DNA hybrids as electrochemical probes

water molecules, and thus missing the direct correlations between water, ONA, and polymer. These solvent based limitations motivate/support development of coarse-grained models for explicit water around biomolecules as well as hybrid or adaptive resolution (atomistic coarse-grained) methods for treating solvent in atomistic/higher detail near the areas of interest within the simulation than elsewhere in the simulation box.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.6b06970. Further details of the all-atom force field parameters, atomistic simulation protocol, metadynamics simulation parameters, and supporting figures for the Results section (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Phone: 302 831 8682. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors thank National Science Foundation grant NSFDMR 1420736 for financial support. This work was also supported in part through the use of information technologies resources at the University of Delaware, specifically the Farber high-performance computing resources.

■

REFERENCES

(1) Watson, J. D.; Crick, F. H. Molecular structure of nucleic acids. Nature 1953, 171 (4356), 737−738. (2) Metzker, M. L. Sequencing technologiesthe next generation. Nat. Rev. Genet. 2010, 11 (1), 31−46. (3) Chmielecki, J.; Meyerson, M. DNA sequencing of cancer: what have we learned? Annu. Rev. Med. 2014, 65, 63−79. (4) Liang, H.; Zhang, X.-B.; Lv, Y.; Gong, L.; Wang, R.; Zhu, X.; Yang, R.; Tan, W. Functional DNA-containing nanomaterials: cellular applications in biosensing, imaging, and targeted therapy. Acc. Chem. Res. 2014, 47 (6), 1891−1901. (5) Tjong, V.; Tang, L.; Zauscher, S.; Chilkoti, A. Smart” DNA interfaces. Chem. Soc. Rev. 2014, 43 (5), 1612−1626. (6) Tay, C. Y.; Yuan, L.; Leong, D. T. Nature-Inspired DNA Nanosensor for Real Time In Situ Detection of mRNA in Living Cells. ACS Nano 2015, 9, 5609. (7) Tang, Y.; Ge, B.; Sen, D.; Yu, H.-Z. Functional DNA switches: rational design and electrochemical signaling. Chem. Soc. Rev. 2014, 43 (2), 518−529. (8) Bandy, T. J.; Brewer, A.; Burns, J. R.; Marth, G.; Nguyen, T.; Stulz, E. DNA as supramolecular scaffold for functional molecules: progress in DNA nanotechnology. Chem. Soc. Rev. 2011, 40 (1), 138− 148. (9) Wang, F.; Lu, C.-H.; Willner, I. From cascaded catalytic nucleic acids to enzyme−DNA nanostructures: controlling reactivity, sensing, logic operations, and assembly of complex structures. Chem. Rev. 2014, 114 (5), 2881−2941. (10) Jones, M. R.; Seeman, N. C.; Mirkin, C. A. Programmable materials and the nature of the DNA bond. Science 2015, 347 (6224), 1260901. (11) Wang, Z. G.; Ding, B. DNA-Based Self-Assembly for Functional Nanomaterials. Adv. Mater. 2013, 25 (28), 3905−3914. (12) Kuzuya, A.; Ohya, Y. Nanomechanical molecular devices made of DNA origami. Acc. Chem. Res. 2014, 47 (6), 1742−1749. J


Article

The Journal of Physical Chemistry B for the detection of DNA. J. Am. Chem. Soc. 2005, 127 (4), 1170− 1178. (34) Jeong, J. H.; Kim, S. W.; Park, T. G. Novel intracellular delivery system of antisense oligonucleotide by self-assembled hybrid micelles composed of DNA/PEG conjugate and cationic fusogenic peptide. Bioconjugate Chem. 2003, 14 (2), 473−479. (35) Rosi, N. L.; Giljohann, D. A.; Thaxton, C. S.; Lytton-Jean, A. K.; Han, M. S.; Mirkin, C. A. Oligonucleotide-modified gold nanoparticles for intracellular gene regulation. Science 2006, 312 (5776), 1027− 1030. (36) Seferos, D. S.; Giljohann, D. A.; Hill, H. D.; Prigodich, A. E.; Mirkin, C. A. Nano-flares: probes for transfection and mRNA detection in living cells. J. Am. Chem. Soc. 2007, 129 (50), 15477− 15479. (37) Hong, B. J.; Eryazici, I.; Bleher, R.; Thaner, R. V.; Mirkin, C. A.; Nguyen, S. T. Directed Assembly of Nucleic Acid-Based Polymeric Nanoparticles from Molecular Tetravalent Cores. J. Am. Chem. Soc. 2015, 137 (25), 8184−8191. (38) Jeong, J. H.; Park, T. G. Novel polymer-DNA hybrid polymeric micelles composed of hydrophobic poly (D, L-lactic-co-glycolic acid) and hydrophilic oligonucleotides. Bioconjugate Chem. 2001, 12 (6), 917−923. (39) Jiang, F. X.; Yurke, B.; Firestein, B. L.; Langrana, N. A. Neurite outgrowth on a DNA crosslinked hydrogel with tunable stiffnesses. Ann. Biomed. Eng. 2008, 36 (9), 1565−1579. (40) Tang, H.; Duan, X.; Feng, X.; Liu, L.; Wang, S.; Li, Y.; Zhu, D. Fluorescent DNA−poly (phenylenevinylene) hybrid hydrogels for monitoring drug release. Chem. Commun. 2009, 6, 641−643. (41) Cao, A.; Tang, Y.; Liu, Y.; Yuan, H.; Liu, L. A strategy for antimicrobial regulation based on fluorescent conjugated oligomer− DNA hybrid hydrogels. Chem. Commun. 2013, 49 (49), 5574−5576. (42) Zhou, L.; Chen, C.; Ren, J.; Qu, X. Towards intelligent bioreactor systems: triggering the release and mixing of compounds based on DNA-functionalized hybrid hydrogel. Chem. Commun. 2014, 50 (71), 10255−10257. (43) Wu, Y.; Li, C.; Boldt, F.; Wang, Y.; Kuan, S. L.; Tran, T. T.; Mikhalevich, V.; Förtsch, C.; Barth, H.; Yang, Z. Programmable protein−DNA hybrid hydrogels for the immobilization and release of functional proteins. Chem. Commun. 2014, 50 (93), 14620−14622. (44) Wan, L.; Chen, Q.; Liu, J.; Yang, X.; Huang, J.; Li, L.; Guo, X.; Zhang, J.; Wang, K. Programmable Self-Assembly of DNA−Protein Hybrid Hydrogel for Enzyme Encapsulation with Enhanced Biological Stability. Biomacromolecules 2016, 17 (4), 1543−1550. (45) Xi, W.; Pattanayak, S.; Wang, C.; Fairbanks, B.; Gong, T.; Wagner, J.; Kloxin, C. J.; Bowman, C. N. Clickable Nucleic Acids: Sequence-Controlled Periodic Copolymer/Oligomer Synthesis by Orthogonal Thiol-X Reactions. Angew. Chem., Int. Ed. 2015, 54 (48), 14462−14467. (46) Mori, T.; Maeda, M. Temperature-responsive formation of colloidal nanoparticles from poly (N-isopropylacrylamide) grafted with single-stranded DNA. Langmuir 2004, 20 (2), 313−319. (47) Alemdaroglu, F. E.; Safak, M.; Wang, J.; Berger, R.; Herrmann, A. DNA multiblock copolymers. Chem. Commun. 2007, 13, 1358− 1359. (48) Seeman, N. C. Nanomaterials based on DNA. Annu. Rev. Biochem. 2010, 79, 65−87. (49) Qi, H.; Ghodousi, M.; Du, Y.; Grun, C.; Bae, H.; Yin, P.; Khademhosseini, A. DNA-directed self-assembly of shape-controlled hydrogels. Nat. Commun. 2013, 4, 2275. (50) Takahashi, K.; Matsuo, M.; Banno, T.; Toyota, T. Micrometersized network structure of novel DNA−lipid conjugates induced by heat stimulation. Soft Matter 2015, 11 (35), 7053−7058. (51) Meng, Y.-F.; Wei, J.; Gao, P.-C.; Jiang, Y. Self-assembling amphiphilic poly (propargyl methacrylate) grafted DNA copolymers into multi-strand helices. Soft Matter 2015, 11 (28), 5610−5613. (52) Prytkova, T. R.; Eryazici, I.; Stepp, B.; Nguyen, S.-B.; Schatz, G. C. DNA Melting in Small-Molecule− DNA-Hybrid Dimer Structures: Experimental Characterization and Coarse-Grained Molecular Dynamics Simulations. J. Phys. Chem. B 2010, 114 (8), 2627−2634.

(53) Jin, R.; Wu, G.; Li, Z.; Mirkin, C. A.; Schatz, G. C. What controls the melting properties of DNA-linked gold nanoparticle assemblies? J. Am. Chem. Soc. 2003, 125 (6), 1643−1654. (54) Park, S. Y.; Stroud, D. Structure formation, melting, and optical properties of gold/DNA nanocomposites: effects of relaxation time. Phys. Rev. B: Condens. Matter Mater. Phys. 2003, 68 (22), 224201. (55) Shchepinov, M. S.; Mir, K. U.; Elder, J. K.; Frank-Kamenetskii, M. D.; Southern, E. M. Oligonucleotide dendrimers: stable nanostructures. Nucleic Acids Res. 1999, 27 (15), 3035−3041. (56) Greschner, A. A.; Toader, V.; Sleiman, H. F. The Role of Organic Linkers in Directing DNA Self-Assembly and Significantly Stabilizing DNA Duplexes. J. Am. Chem. Soc. 2012, 134 (35), 14382− 14389. (57) Gibbs-Davis, J. M.; Schatz, G. C.; Nguyen, S. T. Sharp melting transitions in DNA hybrids without aggregate dissolution: Proof of neighboring-duplex cooperativity. J. Am. Chem. Soc. 2007, 129 (50), 15535−15540. (58) Park, S. Y.; Gibbs-Davis, J. M.; Nguyen, S. T.; Schatz, G. C. Sharp melting in DNA-linked nanostructure systems: thermodynamic models of DNA-linked polymers. J. Phys. Chem. B 2007, 111 (30), 8785−8791. (59) Pechar, M.; Kopečková, P.; Joss, L.; Kopeček, J. Associative diblock copolymers of poly (ethylene glycol) and coiled-coil peptides. Macromol. Biosci. 2002, 2 (5), 199−206. (60) Vandermeulen, G. W.; Tziatzios, C.; Klok, H.-A. Reversible selforganization of poly (ethylene glycol)-based hybrid block copolymers mediated by a de novo four-stranded α-helical coiled coil motif. Macromolecules 2003, 36 (11), 4107−4114. (61) Choi, Y. Y.; Jang, J. H.; Park, M. H.; Choi, B. G.; Chi, B.; Jeong, B. Block length affects secondary structure, nanoassembly and thermosensitivity of poly (ethylene glycol)-poly (l-alanine) block copolymers. J. Mater. Chem. 2010, 20 (17), 3416−3421. (62) Shu, J. Y.; Lund, R.; Xu, T. Solution structural characterization of coiled-coil peptide−polymer side-conjugates. Biomacromolecules 2012, 13 (6), 1945−1955. (63) Zanuy, D.; Hamley, I. W.; Alemán, C. Modeling the tetraphenylalanine-PEG hybrid amphiphile: from DFT calculations on the peptide to molecular dynamics simulations on the conjugate. J. Phys. Chem. B 2011, 115 (28), 8937−8946. (64) Yang, C.; Lu, D.; Liu, Z. How PEGylation enhances the stability and potency of insulin: a molecular dynamics simulation. Biochemistry 2011, 50 (13), 2585−2593. (65) Xue, Y.; O’Mara, M. L.; Surawski, P. P.; Trau, M.; Mark, A. E. Effect of poly (ethylene glycol) (PEG) spacers on the conformational properties of small peptides: a molecular dynamics study. Langmuir 2011, 27 (1), 296−303. (66) Jain, A.; Ashbaugh, H. S. Helix stabilization of poly (ethylene glycol)−peptide conjugates. Biomacromolecules 2011, 12 (7), 2729− 2734. (67) 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 (11), 4053−4060. (68) Woo, S. Y.; Lee, H. Molecular dynamics studies of PEGylated αhelical coiled coils and their self-assembled micelles. Langmuir 2014, 30 (29), 8848−8855. (69) Hamed, E.; Ma, D.; Keten, S. Multiple PEG Chains Attached onto the Surface of a Helix Bundle: Conformations and Implications. ACS Biomater. Sci. Eng. 2015, 1 (2), 79−84. (70) Stanzione, F.; Jayaraman, A. Computational Design of Oligopeptide Containing Poly (ethylene glycol) Brushes for StimuliResponsive Drug Delivery. J. Phys. Chem. B 2015, 119 (42), 13309− 13320. (71) Earl, D. J.; Deem, M. W. Parallel tempering: Theory, applications, and new perspectives. Phys. Chem. Chem. Phys. 2005, 7 (23), 3910−3916. (72) Barducci, A.; Bonomi, M.; Parrinello, M. Metadynamics. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2011, 1 (5), 826−843. K


Article

The Journal of Physical Chemistry B (73) Nielsen, P. E.; Egholm, M.; Berg, R. H.; Buchardt, O. Sequenceselective recognition of DNA by strand displacement with a thyminesubstituted polyamide. Science 1991, 254 (5037), 1497−1500. (74) Lindahl, E.; Hess, B.; Van Der Spoel, D. GROMACS 3.0: a package for molecular simulation and trajectory analysis. J. Mol. Model. 2001, 7 (8), 306−317. (75) Stogryn, A. Equations for calculating the dielectric constant of saline water (correspondence). IEEE Trans. Microwave Theory Tech. 1971, 19 (8), 733−736. (76) Maier, J. A.; Martinez, C.; Kasavajhala, K.; Wickstrom, L.; Hauser, K. E.; Simmerling, C. ff14SB: Improving the accuracy of protein side chain and backbone parameters from ff99SB. J. Chem. Theory Comput. 2015, 11 (8), 3696−3713. (77) Pérez, A.; Marchán, I.; Svozil, D.; Sponer, J.; Cheatham, T. E.; Laughton, C. A.; Orozco, M. Refinement of the AMBER force field for nucleic acids: improving the description of α/γ conformers. Biophys. J. 2007, 92 (11), 3817−3829. (78) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 1983, 79 (2), 926−935. (79) Joung, I. S.; Cheatham, T. E., III Determination of alkali and halide monovalent ion parameters for use in explicitly solvated biomolecular simulations. J. Phys. Chem. B 2008, 112 (30), 9020− 9041. (80) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and testing of a general amber force field. J. Comput. Chem. 2004, 25 (9), 1157−1174. (81) Wang, J.; Wang, W.; Kollman, P. A.; Case, D. A. Automatic atom type and bond type perception in molecular mechanical calculations. J. Mol. Graphics Modell. 2006, 25 (2), 247−260. (82) 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.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M. J.; Heyd, J.; Brothers, E. N.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A. P.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, N. J.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, Ö .; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian 09; Gaussian, Inc.: Wallingford, CT, 2009. (83) Barducci, A.; Bussi, G.; Parrinello, M. Well-tempered metadynamics: A smoothly converging and tunable free-energy method. Phys. Rev. Lett. 2008, 100 (2), 020603. (84) Bonomi, M.; Parrinello, M. Enhanced sampling in the welltempered ensemble. Phys. Rev. Lett. 2010, 104 (19), 190601. (85) Deighan, M.; Bonomi, M.; Pfaendtner, J. Efficient simulation of explicitly solvated proteins in the well-tempered ensemble. J. Chem. Theory Comput. 2012, 8 (7), 2189−2192. (86) Deighan, M.; Pfaendtner, J. Exhaustively sampling peptide adsorption with metadynamics. Langmuir 2013, 29 (25), 7999−8009. (87) Plimpton, S. Fast parallel algorithms for short-range molecular dynamics. J. Comput. Phys. 1995, 117 (1), 1−19. (88) Tuckerman, M.; Berne, B. J.; Martyna, G. J. Reversible multiple time scale molecular dynamics. J. Chem. Phys. 1992, 97 (3), 1990− 2001. (89) Hockney, R. W.; Eastwood, J. W. Computer simulation using particles; Taylor and Francis Group: New York, 1988; p 267. (90) Tribello, G. A.; Bonomi, M.; Branduardi, D.; Camilloni, C.; Bussi, G. PLUMED 2: New feathers for an old bird. Comput. Phys. Commun. 2014, 185 (2), 604−613. (91) Abraham, M. J.; Murtola, T.; Schulz, R.; Pá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. (92) Hanwell, M. D.; Curtis, D. E.; Lonie, D. C.; Vandermeersch, T.; Zurek, E.; Hutchison, G. R. Avogadro: an advanced semantic chemical editor, visualization, and analysis platform. J. Cheminf. 2012, 4 (1), 17. (93) Barducci, A.; Bonomi, M.; Parrinello, M. Linking well-tempered metadynamics simulations with experiments. Biophys. J. 2010, 98 (9), L44−L46. (94) Hinckley, D. M.; Freeman, G. S.; Whitmer, J. K.; de Pablo, J. J. An experimentally-informed coarse-grained 3-site-per-nucleotide model of DNA: Structure, thermodynamics, and dynamics of hybridization. J. Chem. Phys. 2013, 139 (14), 144903. (95) Ouldridge, T. E.; Louis, A. A.; Doye, J. P. Extracting bulk properties of self-assembling systems from small simulations. J. Phys.: Condens. Matter 2010, 22 (10), 104102−104113. (96) Bonomi, M.; Barducci, A.; Parrinello, M. Reconstructing the equilibrium Boltzmann distribution from well-tempered metadynamics. J. Comput. Chem. 2009, 30 (11), 1615−1621. (97) Dwight, Z.; Palais, R.; Wittwer, C. T. uMELT: prediction of high-resolution melting curves and dynamic melting profiles of PCR products in a rich web application. Bioinformatics 2011, 27 (7), 1019− 1020. (98) Stanzione, F.; Jayaraman, A. Hybrid Atomistic and CoarseGrained Molecular Dynamics Simulations of Polyethylene Glycol (PEG) in Explicit Water. J. Phys. Chem. B 2016, 120 (17), 4160−4173.

L


Effects of Polymer Conjugation on Hybridization Thermodynamics of

Recommend Documents