Journey of Poly-Nucleotides through OmpF Porin | The Journal of

Article pubs.acs.org/JPCB

Journey of Poly-Nucleotides through OmpF Porin Hamid Hadi-Alijanvand* and Maryam Rouhani Department of Biological Sciences, Institute for Advanced Studies in Basic Sciences (IASBS), Zanjan, 45137-66731, Iran S Supporting Information *

ABSTRACT: OmpF is an abundant porin in many bacteria which attracts attention as a promising biological nanopore for DNA sequencing. We study the interactions of OmpF with pentameric poly-nucleotides (poly-Ns) in silico. The poly-N molecule is forced to translocate through the lumen of OmpF. Subsequently, the structural and dynamical effects of translocation steps on protein and poly-N molecules are explored in detail. The external loops of OmpF are introduced as the main region for discrimination of poly-Ns based on their organic bases. Structural network analyses of OmpF in the presence or absence of poly-Ns characterize special residues in the structural network of porin. These residues pave the way for engineering OmpF protein. The poly-N-specific pattern of OmpF’s local conductance is detected in the current study. Computing the potential of mean force for translocation steps, we define the energetic barrier ahead of poly-N to move through OmpF’s lumen. We suggest that fast translocation of the examined poly-N molecules through OmpF seems unattainable by small external driving forces. Our computational results suggest some abilities for OmpF porin like OmpF’s potential for being used in poly-N sequencing.

■

Each of these monomers has a β-barrel structure constructed from 16 β strands that are connected with 8 external loops and 8 internal turns.21 One of the external loops called the L3 loop folds back into the β-barrel and interacts with residues in the middle of the barrel to create a constriction zone (CZ). The βbarrel structure forms a pore responsible for diffusion of small hydrophilic molecules.22 By maintaining cytoplasmic solute and metabolite concentrations, OmpF plays an important role in bacterial survival.23 The translocation of many molecules such as colicins, antibiotics, and antimicrobial peptides through OmpF porin was studied by experimental or computational approaches.24−32 Recently, it was claimed that OmpF could identify different single-stranded DNA homopolymers on the basis of the amount of blockage of ionic current and the voltage in which the channel was completely blocked.33 In the current study, using all-atom molecular dynamics simulations, we study the passage of single-stranded poly-G DNA (dG) and RNA (rG) as well as poly-C DNA (dC) homopolymers that consist of five nucleotides through the OmpF channel. We study the effects of this passage on structural and dynamical features of both protein and nucleic acid. We identify the effect of different bases and backbone sugars on the interaction between OmpF and poly-nucleotide (poly-N). Important protein residues are introduced as promising residues for future mutational studies. Finally, after computing potential of mean forces (PMFs) for translocation of poly-Ns, we discuss the possibility of passage of the

INTRODUCTION Having the ability of single molecule identification, nanopores have found different applications in detection of specific molecules1 and sequencing single-stranded nucleic acids.2 Sequencing DNA and RNA molecules is of considerable importance in life sciences and in understanding, diagnosis, and therapy of diseases.3 Nanopore sequencing, that minimizes the need for sample preparation and amplification, is a progressive, accurate, fast, and inexpensive sequencing method.2,4−6 Nanopore sequencing usually involves passage of ions from a membrane-embedded nanopore by applying an external electric field which produces an ionic current. When nucleic acids are added to the system as monomers or polymers, they may block the ionic current differently.7,8 Two main types of nanopores used for sequencing and detection of molecules are synthetic solid-state nanopores, including silicon nitride9−11 and graphene,12 and biological nanopores.2 α-Hemolysin7,13 and Mycobacterium smegmatis porin A (MspA)14 are major biological nanopores used for sequencing RNA or single stranded DNA. In contrast to solidstate nanopores, biological nanopores can be manipulated genetically or chemically and reproduced with significantly small variations in pore diameter.1 In the presence of transmembrane potential, nucleic acids pass through αhemolysin and MspA so fast (one nucleotide per μs)2,15 that nucleotide-induced changes in ionic current cannot be detected in the experimental millisecond time scale. Hence, scientists slow down passage of the nucleotide through the channel by different mechanisms.13,16−18 OmpF is an abundant outer membrane protein of Gramnegative bacteria that consists of three identical monomers.19,20 © XXXX American Chemical Society

Received: January 24, 2015 Revised: April 30, 2015

A

DOI: 10.1021/acs.jpcb.5b00763 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

during cv-SMD. Using the first-order cumulant expansion of Jarzynski’s equality (JE) and performing four independent joint MT and SMD simulations (totally 12 ns) for each poly-N, we computed the PMF of poly-N’s translocation:42,43

mentioned homopolymers through OmpF porin by applying an external driving force on the basis of the total free energy of translocation.

■

THEORETICAL METHODS Molecular Dynamics (MD) Simulations. OmpF porin (PDB ID: 2OMF) was embedded in a POPC bilayer membrane which contained 138−140 lipid molecules (depending on the MD system) in each membrane leaflet. The POPC lipid was selected to close the membrane state to the recently published experiment on poly-N−OmpF interactions.33 The TIP3-solvated cubic box contained protein, bilayer membrane, and pentameric nucleic acid (the extended length of pentamer is similar to the length of OmpF’s lumen) and was ionized with KCl. The ion concentration and total charge were set to 0.15 M and 0, respectively. The MD system was assembled by CHARMM-GUI.34−36 We utilized NAMD 2.937 to perform the MD simulations by using CHARMM27 and CHARMM36 force fields.38,39 Particle-mesh Ewald (PME) full electrostatics was used with a 1 Å grid spacing and periodic boundary conditions. To compute the van der Waals potential, we set the cutoff to 12 Å and started the smoothing step at 10 Å. The temperature was kept constant around 298 K by Langevin dynamics. The Langevin piston Nosé−Hoover method was used to maintain fluctuations of pressure around 1 atm in NPT ensemble simulations. Some preparation steps were performed on simulation boxes before the production runs. Briefly, lipid tails were melted around protein, the protein’s atoms were allowed to partly breathe, hydration of lipid tails was prevented, and finally the system (including water, lipid, ions, protein, and poly-N) was allowed to fully relax. These preparation steps needed a total simulation time of 5 ns. The center of mass (COM) of thte minimized structure of pentameric deoxy ribose with cytosine bases (dC) and guanine bases (dG) and of the RNA version of pentameric dG (rG) was placed at least 20 Å away from the COM of OmpF’s entrance. Poly-Ns were pulled into the pore’s opening from their 5′ ends. The 5′ ends of poly-Ns were blocked by a hydroxyl moiety. Computing PMF Using Joint Metadynamics (MT) and Constant Velocity Steered Molecular Dynamics (cvSMD). The steered molecular dynamics (SMD) approach is used for computing the PMF of poly-Ns’ translocation through the α-hemolysin nanopore.40 To enhance the sampling efficiency and to reduce the amount of dissipative work in the cv-SMD method, a combination of metadynamics (MT) and SMD was utilized similar to the driven metadynamics approach.41 The joint MT and SMD approach was the source of ensembles of structures used for further studies (Supporting Information, Supporting Method). After preparation steps, poly-N was forced through OmpF’s lumen by cv-SMD. The closest base’s nitrogen atom to the ribose or deoxy-ribose moiety of 5′ nucleotide was set as the SMD atom. The SMD atom was constrained harmonically with a force constant of 10 kcal/mol/Å2. The SMD atom was moved with a constant velocity (3 × 10−5 Å/fs) along the z-axis of OmpF’s lumen which was perpendicular to the plane of the membrane. We also restrained poly-N from moving out of the pore cylinder. The poly-N molecule was allowed to deviate at most 20 Å from its equilibrated starting structure, and its endto-end distance was restricted in the range from 10 to 50 Å. A membrane-embedded Cα rim of OmpF was restrained along the z-axis to prevent OmpF from drifting out of the membrane

exp( −PMF/kT ) = ⟨exp(−W /kT )⟩

(1)

Here, T stands for the absolute temperature in Kelvin, k notes the Boltzmann constant, W stands for SMD-derived works, and the brackets note the average over the ensemble reproduced from joint MT and SMD simulations. The details of SMDbased PMF calculation are presented in the Supporting Method (Supporting Information). Electrostatic Potential Surface (EPS). We used the APBSmem software 1.04 to compute OmpF’s EPS in an implicit membrane environment.44 The concentrations of counterions were set to 0.15 M, and the dielectric constants of both the protein and membrane were set to 2. The APBS software 1.3 was used to calculate the volumetric map (threedimensional lattice with a value at each lattice point) of ion charge density.45 The average EPS of protein and poly-N during the translocation process was computed by the PMEPot plugin of VMD. This plugin computes the reciprocal sum of smoothed particle-mesh Ewald that represents the smoothed electrostatic potential.46 This tool is designed to estimate the EPS of all frames of the MD trajectory especially for a large system that contains a protein in an explicit membrane. Calculation of Local Conductance. OmpF forms a waterfilled channel in the membrane. It is possible to find various tunnels with different pore radii and lengths that start from OmpF’s opening and extend toward its lumen. The computed properties of tunnels help us to calculate the local conductance (G) of OmpF on the basis of the Hille equation:47,48 G=

1 ⎡ 1 ⎣⎢ 4g(Rc − r) +

1 4g (Rp − r )

N

+ ∑i = 1

⎤ ΔZi ⎥ gπ (Ri − r )2 ⎦

(2)

Ri stands for the radius of a slice of the tunnel with ΔZi thickness. The radius of the hydration shell of the tunnel is represented by r. The bulk conductivity, cytoplasmic radius, and extracellular radius of the channel are noted by g, Rc, and Rp, respectively. Rc was 6.5 Å and Rp was 5.2 Å in our simulations. For computing the local conductance of OmpF in the presence or absence of poly-Ns, we accepted tunnels with a suitable size of bottleneck (at least 1.4 Å) that had a sufficient length to pass the porin’s lumen. The conductance of OmpF in the absence of poly-Ns was smoothed by percentile filter. The resulting curve was deconvulated to distinct Gaussian peaks. The properties of tunnels were computed by CAVER 3.01 software.49,50 Computing PMF Using Adaptive Biasing Force (ABF). Importance sampling is an approach for calculation of free energy. This approach samples the regions of the reaction coordinate (RC) that are hardly discovered by other methods. Due to the existence of several free-energy barriers, the potential energy surface is not sampled effectively in conventional MD calculations. To guarantee the discovery of new regions in potential energy surface, methods are created that are based on the system’s history. Umbrella sampling, MT, and ABF are examples of importance-sampling history-dependent methods.51 Having information about some topological properties of the potential energy surface is necessary for historydependent simulations. However, the adaptive-bias sampling methods, like ABF and adaptive versions of umbrella sampling, B

DOI: 10.1021/acs.jpcb.5b00763 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B are less dependent on prior knowledge about properties of the potential energy surface. ABF is a method used in both sampling and calculation of PMF along the transition path in MD. The free energy of a process whose progress is described by a suitable RC can be defined as the potential derived from the average of force that works in the RC direction. In the ABF method, the system force in the RC direction is computed instantaneously and the time average of this force is calculated continuously. When sampling reaches a threshold, an external biasing force is exerted along the RC to neutralize the average force. This biasing force opens the way to diffusive motion of the biased atom and therefore to more sampling along the RC. Because in the ABF method the instantaneous force is calculated as the desired physical parameter, this method is considered a thermodynamic integration method. If in an ABF calculation RC is divided into bins with specific and suitable dimensions, the distribution of instantaneous force in each bin will be normal. In the ABF method, the system reaches equilibrium just after finding a good estimation of biasing force.51 However, in pulling methods, the system continuously diverges from the equilibrium state. In the ABF method, the free energy is calculated by forces exerted on the biased atom. If we consider one-dimensional RC (r) for progression of the system, we can describe the derivative of free energy with respect to the RC as follows:51 d F (r ) = dr

∂U (r ) 1 ∂S − ∂r β ∂r

SEABF = (c − b)

(5)

BFE = α(ΔVDW) + (βbound ELC bound − βfree ELCfree) + γ(ΔASA)

(6)

Δ stands for the difference between the average of ensembles’ values in bound and free states of poly-N. VDW and ELC represent NAMD-computed van der Waals and electrostatic pair interaction energies, respectively. In the bound state, the pair interaction energies were computed between poly-N and other non-nucleotide components of the system, whereas for the free state the pair interaction energies were computed between poly-N and water and ions. “α” was set to 0.18, and “β” and “γ” were fitting parameters. To find the adjustable parameters, first we constructed a 100 ps ensemble of equilibrated structures for poly-N in bulk solution and for poly-N trapped in the CZ of OmpF. Then, we fitted the LIEbased BFE to the ABF-derived binding free energy between poly-N and CZ. We used the adjusted parameters, β and γ, to estimate BFE values along the translocation pathway. We were aware of the entropic costs and large perturbations in the system during translocation of poly-N and therefore discussed such results carefully. Normal Mode Analysis. ProDy-1.5.1 was used to perform principal component analysis (PCA) and essential dynamics analysis (EDA).58−61 We performed PCA on Cα atoms of protein to infer dynamical properties of protein in different conditions. The first principal mode of OmpF’s motions was selected by setting the RMSD to 2.0 Å from the prime structure. To estimate changes in the dynamics of OmpF’s pore, we first compressed the MD trajectories by PCA. Then, we analyzed the generated conformations along the first normal mode to calculate the radii of the tunnels’ openings by CAVER 3.01 software. EDA was performed to derive the crosscorrelation matrix of OmpF’s residues along poly-N’s translocation path. Structural Features of Poly-N. The X3dna-v2.1 software was used to derive the structural properties of poly-N molecules along the translocation pathway.62,63 These properties that include base stacking, rise, helix twist, twist, roll, inclination, and the set of angles α, β, γ, δ, ε, ζ, and χ were used for clustering. The shared area between adjacent bases was considered as a marker of base stacking. Protein Structural Network (PSN). The PSN approach considers protein as a graph in which amino acids are the graph’s nodes connected via nonbonding interactions. The stability of the structural network is defined by the strength of nonbonding interactions. The strength of an interaction

Here, the first term on the right indicates the physical forces exerted on the system that result from the potential energy function of the system. The second term represents the geometric entropy. β stands for the kT term. In the ABF algorithm, to obtain an estimation of dF(r)/dr, the RC is divided into small bins with specific dimensions. In this equation, dF(r) is an early estimate of free energy as a function of the RC. During simulation, biasing force or fABF is exerted on the biased atom and expressed as fABF = dF(r)/dr. fABF is the negative of the average instantaneous force exerted (⟨f r⟩r) on the biased particle in the direction of the RC:40 d F (r ) = −⟨fr ⟩r dr

(1 + 2κ )

σ2 notes the variance of the instantaneous force, N stands for the total number of force values that are used to construct the ABF-derived PMF, and κ is the correlation length for the series of calculated instantaneous forces.53,54 Linear Interaction Energy (LIE)-Based Binding Free Energy (BFE) Calculation. If we suppose that a system is oscillating harmonically, we can assume a linear response to perturbations. The LIE-based BFE calculation follows such an assumption. It is a simple model of host−guest interactions that uses the ligand’s interaction energies in free and bound states to compute the BFE. An improved version of the LIE-based BFE was used to resolve the computed structural clusters of poly-Ns along poly-N’s translocation path:55−57

(3)

r

σ N

(4) ABF

To avoid the nonequilibrium state after exerting f in the NAMD-implemented ABF algorithm, a good estimate of fABF is found by allowing the system to reach an equilibrium state and to be sampled sufficiently. Then, the calculated fABF is exerted on the system by a ramp function. In this manner, the biased particle gains the self-diffusion property and can even step backward. We split the joint MT and SMD-proposed poly-N translocation pathway into nonoverlapping bins of length 0.15 Å. The PMF was computed along the z-axis of the channel as a function of the distance between the COMs of OmpF and poly-N. To find a reasonable simulation time before applying the biasing force, the ABF code was run in each bin for 0.01, 1, and 5 ps. Eventually, the biasing forces were applied after a 5 ps sampling in each bin. In ABF runs, all of OmpF’s atoms were free but poly-N’s atoms had SMD-like constraints. Each ABF computation was continued for 8 ns and repeated three times in independent simulations. The maximum of standard error between point b and c that resides on the reaction coordinate in ABF-derived PMF (SEABF) is computed by52 C

DOI: 10.1021/acs.jpcb.5b00763 J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

■

Article

RESULTS AND DISCUSSION To facilitate following up the results, we present the results and discussion in four sections: (1) Translocation of poly-Ns through OmpF: This section covers the results of joint MT and SMD simulations. The consequences of translocation of poly-Ns through OmpF’s lumen and the effects of poly-Ns on the channel’s local conductance are discussed. (2) Protein side of Poly-Ns’ translocation process: In this section, focus is on the effects of poly-Ns’ translocation on structural aspects of protein. (3) Poly-N side of poly-Ns’ translocation process: Poly-N experiences many nucleotide-specific structural changes during traveling along the lumen of OmpF. (4) Energetics of the translocation process: On the basis of the energetics of translocation, the possibility of in vitro poly-Ns’ translocation is discussed. Translocation of Poly-Ns through OmpF. Using nanopore structures is one of the revolutionary methods for sequencing poly-Ns. Some channel proteins are suitable candidates for being used as a Coulter counter in the nanoscale. To detect the poly-N-derived blockage of ion current in single channel current recording methods in wet or computational experiments, researchers use an external driving force. The external potential forces poly-N to move through the lumen and therefore partially or completely block the ion current through the channel.71 Analyzing the ion current and conductance of the channel provides a fine tool to determine the sequence of translocating poly-N. As in single channel current recording methods, in the current study, poly-Ns are forced along the nanopore in a nonequilibrium manner by the SMD approach.40 Briefly, a dummy atom is connected to the 5′ end of poly-N via an imaginary spring. An external force is imposed on the spring via the dummy atom to pull poly-N with a constant velocity by cv-SMD.72 An animation presented in the Supporting Information shows poly-N’s translocation through OmpF using joint MT and SMD (Supporting Information, animation). To detect local barriers ahead of poly-N, the amount of local work is derived from the movement of poly-N every 200 fs for dG, rG, and dC poly-N molecules as a function of the SMD atom’s distance from the COM of the protein (Figure 1). The variations in the amount of local work necessary to pull poly-Ns through the channel signify the existence of different barriers in the moving pathway of each poly-N. The required local work for translocation of poly-Ns indicates that energy barriers are distributed in-homogeneously along poly-Ns’ moving pathways. Poly-Ns experience large barriers, which need a high amount of local work to be overcome at the entrance of the constriction zone (CZ). By considering the error bars of local work values, it appears that the amount of local work required to exit from the CZ seems higher for deoxy types of poly-Ns than for rG. It seems that exiting the CZ is a backbone-dependent step. These observations suggest that the type of backbone leads to poly-Ns’ discrimination when they are exiting the CZ. Poly-N molecules are highly hydrated in bulk solution. This hydration plays a critical role in 3D structural polymorphisms of poly-Ns. During translocation of poly-N through the lumen of OmpF, the amount of poly-N’s hydration varies in a nucleotide-dependent manner (Figure 2). There are some hydration-modulator niches (highly populated regions in the

between two amino acids is directly proportional to the number of their atom pairs located within a 4.5 Å distance. There are cutoff values for the strength of interactions.64 A higher cutoff selects stronger interactions between nodes. An interaction with a strength higher than a specified cutoff value is considered a contact. A contact that persisted for at least 25% of the simulation time was assigned as a stable contact. A residue could be a hub if it had at least four stable contacts. If a residue remained a hub in interaction cutoffs 3 to 5, it was counted as a master residue (MR). PSN analysis was performed by Wordom 0.22.65,66 Wordom was also used to count molecules in a 3 Å thick shell around poly-N and to compute poly-N’s solvent accessible surface area along trajectories. The hydration ratio of poly-N was defined as the ratio of the number of water residues to nonwater residues in a 3 Å shell around poly-N. Affinity Propagation Clustering (APC). In clustering methods, usually the input is the amounts of similarity between data points. In k-centers clustering methods, first the clusters’ exemplars are chosen randomly and then they are refined. However, in affinity propagation, all data points are assumed as possible exemplars and the number of clusters is not dictated by the algorithm. In this method, each data point is considered as a node in a network. On the basis of the similarity of data points, messages are exchanged between these nodes and finally stable exemplars and their relative clusters are recognized. These messages include “responsibility” and “availability” of nodes. Responsibility is sent from a data point to a candidate exemplar and indicates the interest of the data point in this candidate to be an exemplar, whereas availability is sent from a possible candidate exemplar to a data point and indicates how much the candidate can be an exemplar for the data point. During clustering, both exemplar and related data point roles are possible for all data points. The responsibility and availability messages are updated during the APC run. Actually, these messages indicate the tendency of a data point to select another data point as its exemplar at each moment. When data points reach an agreement on their exemplars, the algorithm will finish.67 Before performing the APC method, it is necessary to define poly-Ns’ appropriate descriptors. X3dna was used to compute structural metrics of poly-Ns. To reduce the dimension of X3dna-computed structural metrics, the root-mean-square distances of nucleotides’ metrics from the structural metrics of standard DNA polymorph structuresA, B, C, D, E, Z(C), and Z(G)were determined for poly-Ns’ structures during the translocation process.68 The average distance of five constituent nucleotides of poly-N from each standard polymorph structure of DNA was fed into APC 1.1.1 which was implemented in R 2.14 software.69 Applying a Constant Electric Field to the System. To evaluate the ease with which poly-Ns translocate through OmpF, using the NAMD external electric field feature, we forced the dC molecule through porin’s lumen by applying a constant electric field (2 V) to the simulation box perpendicular to the membrane plane. Miscellaneous. The fitting multi Guassian peaks feature of the OriginPro software was used to deconvolute the conductance of poly-N-free OmpF to distinct peaks. The OriginPro software was also used to smooth the data. All molecular images were generated by VMD 1.9.1.70 Simulations were run on three Corei5 PCs. D

DOI: 10.1021/acs.jpcb.5b00763 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 1. Local work for poly-Ns’ translocation and schematic representation of OmpF’s proposed regions. The upper panel indicates the amounts of local work required for moving poly-Ns through OmpF’s lumen as a function of the distance of the SMD atom from the protein’s COM in angstroms. The local work values are smoothed by averaging over the moving window that is 50 ps in width. The gray error bars represent the standard error of mean (SEM) of four independent simulations. In this study, the structure of OmpF is divided into three regions (outer, CZ, and inner) along the z-axis which is perpendicular to the plane of the membrane (lower panel). The positions of the regions’ boundaries are represented in relation to the COM of protein in angstroms. The L3 loop resides in CZ.

Figure 2. Fraction of water molecules in poly-N surrounding residues depicted as a kernel density map. The ratio of water residues to residues that surround poly-N is represented as a function of the distance of the SMD atom from protein’s COM in angstroms. Color bars indicate the population density.

density map) along OmpF’s lumen that dehydrate or hydrate poly-Ns differentially. It seems that poly-N should be dehydrated for entering CZ. dC undergoes an early dehydration before entering CZ. This observation points to the possibility that the early dehydration of poly-N is determined by the base type. It was reported that the hydration of monomeric guanine was more favorable than that of monomeric cytosine.73 Therefore, our observation seems in agreement with the mentioned study. The dehydration phenomenon is usually an indicator of hydrophobic interactions. As dC becomes dehydrated sooner and higher than other poly-Ns do, dC might interact with OmpF or with another part of its own structure by hydrophobic interactions. Poly-N’s translocation is not a single-step process. Possibly, three distinct steps are conceivable for this phenomenon: the incoming molecule is captured by the channel’s entrance, the captured molecule starts the penetration step, and the penetrated molecule passes the lumen in the translocation step. The highest amount of poly-Ns’ dehydration occurs when poly-N passes the CZ during the penetration step.

The conformation and dynamics of the channel’s mouth play a critical role in capturing and penetration steps. The capture rate was a determinative factor in nucleotide differentiation by a protein channel.74,75 Suitable direction of poly-N relative to and sufficient interactions of poly-N with pore’s opening determine successful capturing and penetration steps.76 The potential energy of poly-N’s interactions with the outer region of OmpF and other non-nucleotide atoms in that region is represented as a function of poly-N’s distance from the COM of the protein in Figure 3. The highest interaction energy belongs to poly-N’s interaction with the entrance of the channel. rG and dG have higher interactions in the channel’s mouth than dC does. It seems that base type rather than backbone type plays the determinative role in the interaction energy of the capturing step. The mentioned difference in the interaction energy of polyNs with other parts of the simulation system continues in the entrance of the CZ. The type of poly-Ns’ backbone determines the amount of interaction energy of poly-Ns with other parts of the system during their movement in the CZ. The gradual divergence of rG from dG in the CZ supports the determinative E

DOI: 10.1021/acs.jpcb.5b00763 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

in the lumen would cause higher amounts of ions in OmpF’s lumen than that of dG and we would not observe different ion concentrations in rG- and dG-occupied channels. However, our observations indicate that dC causes the lowest amount of ions in OmpF’s lumen during translocation and surprisingly rG leads to a little change in ion concentration in comparison to poly-Nfree OmpF. Such discrepancies might be rooted in the special 3D structure of poly-N in OmpF’s lumen. There are fine MD simulation reports about translocation of molecules through the lumen of protein channels and the effect of these translocations on channels’ local conductance.46,48 Variations in conductance provide a metric for partial discrimination of molecules that pass the channel. In a convenient model, the fluctuations of local conductance are related to local dynamical changes of the lumen’s tunnels by the Hille equation (see eq 2 in the methods).47 All possible tunnels that traverse the OmpF’s lumen are determined. The length and bottleneck radius of each tunnel is computed. Then, using the Hille formula and the tunnels’ information, we estimate the local conductance of OmpF in the presence or absence of poly-Ns. The smoothed values of local conductance are represented as a function of poly-N’s distance from the COM of porin (Figure 5).

Figure 3. Definition of the interaction energy of poly-N with the barriers ahead. The pair interaction energy of poly-N with nonnucleotide parts of the MD system (including water molecules, ions, and the atoms of the membrane and protein before CZ) is represented as a function of the SMD atom’s distance in angstroms from the protein’s COM. The gray error bars represent the SEM of four independent simulations. The vertical axis values are smoothed by the adjacent-averaging method.

role of backbone type in the interaction energy of poly-N with protein within the CZ. Differences in the interaction energy of poly-Ns with other parts of the system, variations in required local work for translocation, and specific amounts of dehydration during translocation of poly-Ns suggest that the OmpF channel may check poly-N’s base type at the pre-CZ region; the type of poly-N backbone is verified later after reaching the CZ region. Perturbation of the channel’s ionic current is the main signal to distinguish base types of poly-Ns in the Coulter counter approach. The ion concentration in the lumen may be an indicator of the effect of poly-N’s translocation on the ionic current. The distribution of ion concentrations in the lumen of OmpF along the simulation time is depicted in Figure 4. The

Figure 5. Variations in local conductance of OmpF in the presence of poly-Ns. The local conductance of OmpF in the presence or absence of poly-N is computed using eq 2. The relative conductance is calculated by dividing the local conductance by the maximum local conductance of poly-N-free OmpF. The relative conductance of polyN-free OmpF is depicted in the inset as a function of simulation time. Porin’s conductance is fitted to five distinct Guassian peaks named P1−P5. Gray error bars indicate the SEM of three independent 5 ns simulations. In the main figure, the relative conductance of OmpF in the presence of each poly-N is represented as a function of the distance of the SMD atom from the COM of OmpF in angstroms. The gray error bars represent the SEM of four independent simulations. The vertical axis values are smoothed by setting the percentile filter to 85%.

Figure 4. Concentrations of ions in the lumen of OmpF. The histograms of ion concentrations are computed for ensembles. The numbers of ions and water molecules that reside in the lumen are counted every 100 ps, and the data are used to construct the histograms. The gray error bars represent the SEM of four independent simulations.

We compute the maximal conductance of OmpF monomer in poly-N-free MD simulations to be about 220 pS which is close to experimental measurements (313 or 350 pS) at the same ion concentrations.77,78 Analysis of OmpF’s conductance indicates innate ordered fluctuations in porin’s conductance (Figure 5, inset). The frequency of innate beats is about 1/ns. The observed beats may originate from the cross-talk between OmpF porin and the lipid environment at the molecular level.79 It was proposed that salt concentration is able to modulate OmpF current.80 The local concentration of ions in the

reference data is the distribution of ion concentrations of the lumen in the nucleotide-free simulations of OmpF porin. The largest decrement in the concentration of the lumen’s ions due to translocation of poly-N belongs to dC, whereas rG has the lowest effect on ion concentrations. If the dimensions of poly-N’s bases were the only decisive factor for the amount of ions in the lumen, the presence of dC F

DOI: 10.1021/acs.jpcb.5b00763 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 6. Distribution of tunnels’ radii in the pore region of OmpF. After computing principal modes of joint MT and SMD trajectories for the outer region of OmpF, the essential modes of the region are extracted. Ensembles of structures are obtained and analyzed to find possible tunnels that start from the outer region. The bottleneck radii of computed tunnels are measured, and the distribution of radii in angstroms is presented. The gray error bars represent the SEM of three or four independent simulations in the absence or presence of poly-N, respectively.

more than dG did in an experiment.33 This observation may hint at our result. Protein Side of Poly-Ns’ Translocation Process. The interactions of poly-Ns with various parts of OmpF induce many structural and dynamical changes in protein which will be illustrated in the current section. If we assume that the native state of a protein is an ensemble of dynamically interconverting structures around a central structure which is in global energy minima, it is possible to consider a more realistic definition of the native state, that is, an ensemble of native structures. Therefore, just a thermodynamic view is not sufficient to illustrate the state of molecules and a dynamical perspective is necessary and acts in a complementary fashion. Purine-containing poly-Ns interact tightly with porin’s opening (Figure 3). These interactions trigger the capturing step which causes many changes in the dynamics of the pore region. Analysis of the MD trajectory-derived motions of OmpF by principal component analysis (PCA) and essential dynamics analysis (EDA) alleviates the need for the harmonic motion assumption in the elastic network model of protein dynamics. An ensemble of structures defined by the first (slowest) principal mode is constructed for the pore region of OmpF by analyzing joint MT and SMD trajectories. To represent the dynamics of the pore region, the radii of tunnels starting from the openings of the collective mode-derived OmpF structures are measured and their distribution is computed in the presence and absence of poly-Ns (Figure 6). The distribution of tunnels’ radii indicates that each poly-N induces changes in the dynamics of the pore region in a different manner. dC restricts the distribution of tunnels’ radii to lower values with respect to poly-N-free OmpF. Guaninecontaining poly-Ns induce a population of tunnels with radii higher than the tunnels’ radii in the absence of poly-N. Simultaneous attention to Figures 5 and 6 reveals that the higher local conductance of dG-treated OmpF in the pre-CZ region than that of poly-N-free OmpF might be a consequence of induction of tunnels with wide openings. However, dC induces many tunnels with a smaller size of pore in comparison to free OmpF which may retain a small amount of local conductance when dC is in the CZ. Such narrow tunnels may create leaky paths and allow a shallow conductance.

membrane−protein interface may be the source of innate beats in poly-N-free OmpF. The frequency and amplitude of innate beats of OmpF are possibly affected by environmental variables such as membrane potential or interaction with a ligand. Studying the local conductance of the OmpF channel as a function of poly-N’s distance from the COM of protein reveals that poly-N molecules change the local conductance of the channel in a specific manner. As dG moves toward OmpF’s opening, the conductance increases compared to the conductance of the poly-N-free channel. During penetration of dG, the conductance decreases continuously (Figure 5). Movement of rG and dC toward the CZ results in a mild decrease of conductance in comparison to free OmpF. On the basis of these observations, poly-Ns affect OmpF’s local conductance during traveling toward the CZ. These observations indicate that poly-N is decreasing porin’s conductance before reaching the CZ. Therefore, it seems that the conductance-based discrimination among poly-Ns is started when poly-N interacts with the pre-CZ region of OmpF. The three types of poly-N in the current study suppress OmpF’s local conductance differentially in the CZ. These results suggest that the time of suppressing the local conductance of OmpF to the minimum might be an indicator of the type of poly-N; dG causes the minimum current sooner than rG does, and rG shows this effect sooner than dC does (Figure 5). The dynamics of poly-Ns when they are in close contact with the L3 loop in the CZ might lead to a residual conductance. The remaining conductance of OmpF after traveling of dC through the CZ is smaller than the remaining conductance after translocation of dG and rG. When poly-Ns leave the inner region of OmpF, the porin’s local conductance starts recovery to poly-N-free OmpF’s conductance in a poly-N-dependent manner. The rate of local conductance recovery is as follows: rG > dC > dG. The presence of poly-N affects the distribution of local conductance values extracted from MD simulations (Supporting Information, Figure S1). The distribution of local conductance values indicates that the translocation of dC through OmpF leads to maximum suppression of porin’s conductance. The dG molecule decreases OmpF’s conductance less than dC and more than rG. dC blocked the OmpF current G

DOI: 10.1021/acs.jpcb.5b00763 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

the proposed leakiness in OmpF’s local conductance when dC is present at the CZ (Figure 5). By essential dynamics analysis (EDA) of OmpF’s MD trajectories in the presence and absence of poly-N, the crosscorrelation heat maps of residues’ fluctuations are calculated (Supporting Information, Figure S2). In EDA, the crosscorrelation matrix is the residues’ fluctuation covariance matrix which considers the amplitude of changes in dynamics and discards the direction of changes in fluctuations of residues. The study of cross-correlation matrices points to the fact that distinct regions on cross-correlation heat maps appear or fade by translocation of each type of poly-N through OmpF’s lumen. Such a study guides us to explore the changes in OmpF’s structural network due to poly-Ns’ translocation. The structural network approach assumes that the protein structure network (PSN) is a network of amino acid nodes which are connected via nonbonding interaction links. The nodes which have strong links with surrounding nodes and are essential for structural networks’ stability create hubs.64 The strength of contacts is determined by the interaction cutoff parameter; higher interaction cutoff selects stronger noncovalent contacts among amino acids in a network. Properties of hubs are altered by specific interactions that appear or disappear in OmpF due to the interaction of poly-N with distinct regions of protein. We find various hubs for OmpF in the presence or absence of poly-Ns in different interaction cutoffs. We define condition-specific master residues (MRs) by selecting the residues that remain hubs in cutoffs 3 to 5. Some master residues appear by the interaction of any poly-Ns with OmpF (nucleotide MRs), some appear in response to the interaction of a specific type of poly-N with protein (poly-N-specific MRs), and some MRs become visible only in a poly-N-free porin (protein MRs) (Table 2). These results reveal that OmpF

The interactions of poly-N molecules with porin do not change the dynamics of OmpF’s pore exclusively, and the presence of such nucleic acid molecules induces serious dynamical changes in all parts of OmpF. After performing PCA for all Cα atoms of protein in trajectories, the first slowest (nonzero) mode of the OmpF structure is derived in the presence or absence of poly-N molecules. Such calculations determine how translocation of poly-N through porin may affect the dynamics of OmpF. We summarize the changes in dynamics of residues in the presence of poly-Ns at distinct regions of porin’s structure in comparison to the nucleotidefree structure of OmpF in Table 1. If the changes in the Table 1. Effects of Poly-N’s Translocation on the Dynamics of Different Sections of OmpF outer

CZ

inner

↑a ↓ ≈ ↑ ↓ ≈ ↑ ↓ ≈

dG

dC

rG

0.23b 0.40 0.36 0.34 0.32 0.34 0.48 0.31 0.21

0.34 0.33 0.33 0.42 0.28 0.30 0.58 0.10 0.32

0.16 0.60 0.24 0.17 0.46 0.37 0.24 0.31 0.45

a

The changes in the dynamics of amino acids are computed by comparing PCA values (fluctuations) of amino acids in the presence of poly-Ns with corresponding PCA values in the absence of poly-Ns. The upward arrow stands for an increment of at least 50% in PCA value in the presence of poly-N compared to the corresponding value in the absence of poly-N. The downward arrow indicates a decrement of at least 50% in PCA values. The ≈ symbol demonstrates changes that do not reach the 50% cutoff and means no significant dynamical change compared to poly-N-free OmpF. bThe values are the fraction of amino acids with corresponding behavior (increment, decrement, or no change in dynamics in response to the presence of poly-Ns) in each section of OmpF.

Table 2. MRs of OmpF That Are Derived from PSN Analysisd

dynamics of a residue are greater than 50% of its dynamics in nucleotide-free porin, the residue is counted as a residue whose dynamics is affected by the presence of poly-N. PCA of protein motions indicates that dG and rG are able to suppress the dynamics of at least 40% of residues of the OmpF’s outer region (Table 1). The pattern of changes in residues’ dynamics becomes different at the CZ region of OmpF. The dC polymer induces high fluctuations in CZ residues, whereas rG suppresses the dynamics of most of the CZ residues. After passing the outer region and the CZ, poly-N reaches the inner region of OmpF. A poly-N at the inner region also affects the dynamics of residues in a specific manner: dG and dC poly-Ns induce fluctuations in at least 48% of residues of the inner region; however, rG does not change the dynamics of most of these residues. These results imply that the base type of poly-N defines the changes of dynamics in the outer region of OmpF, whereas the backbone type of poly-N affects the dynamics of the CZ and inner regions of OmpF. The reduction of protein’s dynamics may partially be rooted in higher interaction of poly-Ns with different regions of OmpF. Guanine-containing poly-Ns have high interactions with the outer region of OmpF (Figure 3), and such interactions may reduce the dynamics of the pore. On the other hand, the interaction of dC with OmpF enhances the dynamics of the CZ region’s residues. Such an increase in dynamics may rationalize

a

Protein MRs are special residues that persist in the PSN of OmpF in the absence of poly-Ns. bPoly-N-specific MRs are residues that remain in the PSN of porin in response to a specific type of poly-N and do not appear in other conditions. cNucleotide MRs are special master residues induced in OmpF in response to all types of poly-Ns. Residues that reside in the outer, CZ, or inner regions of OmpF are colored orange, blue, or green, respectively. dThe master residues of OmpF and their positions in OmpF are defined.

responds to the presence of different types of poly-Ns by changing its structural networks. The provided MRs might be suitable choices to modulate the response of OmpF to specific types of poly-Ns in future mutational studies. The dC poly-N induces widespread changes in the PSN and creates specific structural networks at the outer, CZ, and inner regions of OmpF (Table 2). dG stimulates the appearance of H

DOI: 10.1021/acs.jpcb.5b00763 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 7. Affinity propagation clustering of poly-N structures during the translocation process. Poly-N structures are segregated into distinct clusters by the APC approach. The APC results of dC, dG, and rG are represented in panels a, b, and c, respectively. The exemplar species of each cluster resides in the center of the cluster’s graph. In the common form of representing APC results, the distance of a cluster’s member from its corresponding exemplar is represented by the length of the member−exemplar interconnecting line. Clusters are colored differentially. The position I

DOI: 10.1021/acs.jpcb.5b00763 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B Figure 7. continued

of each cluster along the translocation pathway is defined on the horizontal axis. The horizontal axis denotes the distance of the cluster’s member from the COM of OmpF in angstroms. The vertical axis represents the LIE-based BFE of the cluster’s member in kcal/mol. The structural parameters of exemplars are summarized in microarray-like tables adjacent to each cluster. The exemplar number and its position of appearance along the translocation pathway are shown at the top and bottom of each table, respectively. The columns of the tables from left to right represent the nucleotides of each exemplar from 5′ to 3′. The 1st to 11th rows demonstrate helix-twist, roll, twist, inclination, α, β, γ, δ, ε, ζ, and χ parameters, respectively. Therefore, each table characterizes the structure of an exemplar. The color codes of the structural metrics values are imbedded at the bottom of the figure.

MRs at the outer region. rG induces MRs in the outer region and CZ. Therefore, OmpF constructs special MRs in the outer region in response to the presence of poly-Ns. Such specific MRs change protein dynamics in response to translocation of poly-Ns and may regulate fluctuations of the pore’s size in a poly-N-specific manner (Figure 6). Focusing on protein MRs indicates that more than 50% of these MRs reside in the L3 loop. This observation implies that the L3 loop plays a critical role in poly-N-independent dynamics of OmpF. Other vital MRs in OmpF are nucleotide MRs. The nucleotide MRs, Y32 and F118, are independent of the type of poly-N, and these residues persist as MRs whenever any of the poly-Ns pass through the lumen of OmpF. The mentioned aromatic residues act as a hydrophobic snap fastener that fastens the L3 loop to the wall of OmpF and restricts the motion of the L3 loop. The pair interaction energy between Y32 and F118 is computed in both poly-N-free and poly-Ncontaining OmpF. In the absence of poly-Ns, such a snap is loose and therefore L3 moves freely (Supporting Information, Figure S3). However, when a poly-N reaches porin’s opening, such a snap becomes tight and restricts the motion of the L3 loop by fixing the loop to the wall of OmpF. Poly-N Side of Poly-Ns’ Translocation Process. Because OmpF’s lumen seems like a confined environment with physical properties different from that of the bulk solvent milieu,81,82 studying the dynamical and structural properties of poly-Ns within OmpF’s lumen is not as straightforward as exploring the dynamics of poly-N at bulk solvent. Therefore, we use a clustering approach to figure out possible structurally orchestrated populations of each poly-N while it is moving through the lumen of OmpF. Clustering procedures need comprehensive descriptors of the target population, so we compute many base-dependent parameters of nucleic acids, such as rolling, twist, and stacking, and many backbonedependent metrics, like gamma and epsilon angles and helixtwist. The affinity propagation clustering (APC) method was used successfully in various fields such as analysis of the trajectories of protein folding.83,84 After computing appropriate descriptors for poly-N structures during translocation, we perform APC. The resulting clusters and clusters’ representatives are determined and depicted in the panels of Figure 7. The results of structure-based clustering are represented as a function of poly-N’s distance from the protein’s center of mass and the binding free energy (BFE) between poly-N and OmpF. To clearly resolve clusters’ positions, the linear interaction energy (LIE)-based BFE is computed for each poly-N along the translocation pathway (see eq 6 in the methods). Attention to BFE of members in a cluster indicates that there is a spectrum of BFE values in each cluster. Because the members of a cluster are gathered on the basis of structural similarities, the spectrum of BFE values in each cluster might imply that the protein side of poly-N−protein interaction has considerable plasticity at the interaction site. The existence of such distinct clusters implies

the existence of specific types of poly-N−protein interactions, the spatial effects of lumen’s microenvironments on the traversing poly-N, and the effects of resident ions and water molecules of the lumen on the dynamics of poly-N. dC and dG polymers tend to have two different populations of structural species at the CZ, whereas rG has a single population there. Deoxy poly-Ns segregate into two different populations, whereas rG shows three distinct structural species before reaching the CZ. The results of APC indicate that contents of OmpF’s lumen prepare the way for segregation of poly-N structures to individual populations along the translocation path. About the rG molecule, the affinity between the representative structure of the cluster and OmpF is higher at the CZ than at the outer and inner regions. BFEs show a decreasing pattern from the outer region to the inner section of OmpF for representative structures of dG. The BFEs of representative structures of dC show an inverse trend. These observations imply that protein−poly-N interaction becomes tighter for dG and looser for dC by movement of poly-N from the outer to inner regions of the pore. Changes in structural parameters of the representative structure of each cluster are summarized in separate tables beside each cluster in Figure 7. Changes in the metrics used for clustering poly-N structures are divided into two distinct categories: (1) the changes in parameters related to adjacent bases (interbase) which are represented in rows 1 to 4 of tables of clusters and (2) the changes of parameters that describe the characteristics of discrete bases (intrabase) which are noted in rows 5−11 of tables of clusters. The most prominent variation in the metrics of rG occurs for the interbase parameters of the CZ-resided cluster in comparison to the adjacent clusters, whereas for dG this distinct variation occurs for the intrabase parameters of CZ-resided clusters. This dissimilarity indicates that the ribose of poly-N’s backbone induces conversion of intrabase changes to interbase changes of parameters. Such a clear pattern of interbase or intrabase changes of parameters is not obvious for the CZ-resided clusters of dC. It seems that variations in intrabase parameters stimulate changes in base stacking at the CZ. For example, the representative structure of the CZ-resided clusters of dG shows prominent changes in intrabase parameters and so a weak stacking there in comparison to adjacent sections (Supporting Information, Figure S4). However, in dC and rG whose changes in interbase parameters increase respectively, the amount of stacking at the CZ is high and increases accordingly. There is an obvious difference between the base stacking of dG and that of rG when these poly-Ns are traversing the CZ (Figure S4, Supporting Information). It seems that OmpF facilitates the persistence of base stacking of rG in the CZ. The deoxy ribose-containing poly-Ns lose their base stacking at the CZ in part. These observations again point to the idea that the CZ’s environment is able to discriminate between ribose and J

DOI: 10.1021/acs.jpcb.5b00763 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

entities that connect nonequilibrium work (derived from varieties of SMD methods) to the equilibrium-based free energy (see eq 1 in the methods).42,88−91 Using the JE, we obtain a primary estimation of the PMF of poly-N’s translocation from joint MT and SMD simulations (Figure 8).

deoxy ribose poly-N molecules; it is possible that CZ checks the backbone of poly-N. In spite of similar base dimensions in dG and rG, these poly-Ns do not show the same base stacking at the CZ. It seems that the OH group of rG’s backbone has a critical role in base stacking especially at the CZ. We expect higher stacking between large bases of dG than between small bases of dC. However, backbone and base rotations lead to higher stacking of dC than that of dG at the CZ. The snapshots of OmpF’s tunnels and poly-Ns’ structures in the CZ during translocation are presented in the Supporting Information, Figure S5. The mentioned results describe structural and dynamical changes in poly-N and OmpF during translocation of poly-N through OmpF. Those parts show the distinct effects of mutual OmpF−poly-N interactions but do not discuss the easiness of poly-N’s translocation through OmpF. Energetics of the Translocation Process. The computation of the electrostatic potential surface (EPS) is a primary tool to detect ligand−protein interactions. Complementary EPSs of protein and ligand suggest their possible interaction. For example, many DNA-binding proteins have a positive EPS in their DNA-binding surface which complements the negative EPS of DNA.85−87 To define the EPS of OmpF, we need to perform Poisson−Boltzmann electrostatics calculations in the membrane milieu by APBSmem.44 APBSmem calculates the EPS of proteins in an implicit membrane. The EPS of OmpF is computed in the presence and absence of poly-Ns at the CZ (Supporting Information, Figure S6). These results indicate that, in spite of extensive negative EPS in the membranous part of OmpF, there is a narrow positive EPS from the porin’s opening toward the lumen’s wall opposite the L3 loop. Such positive EPS may provide the binding region for molecules with negative EPS such as poly-Ns. The computed densities of ionic charges on OmpF in poly-N-bound and poly-N-free modes reveal a cloud of negative ions that reside on a part of porin’s opening and the pre-CZ region (Figure S6, Supporting Information). These calculations imply that poly-Ns might bind favorably to OmpF from opening to the pre-CZ region. We computed the particle-mesh Ewald (PME)-based EPSs of OmpF and OmpF−poly-N systems using the information on MD simulations. The PME-based EPSs indicate that there is an extensive negative EPS on almost all regions of OmpF and limited positive EPSs at the outer and inner regions of porin. Such a large negative EPS of OmpF may hinder ligands with negative EPS such as poly-Ns (Figure S6, Supporting Information). Therefore, these results weaken the possibility of translocation of poly-N through OmpF. We need a precise metric to study the possibility of translocation of poly-N through OmpF’s lumen such as the potential of mean force (PMF). To find good estimates of the proper external driving force for passing a special type of poly-N through OmpF’s lumen, we need to know the energy barriers of translocation. Computing the PMF of translocation steps partly resolves this issue. Slow SMD simulations are common to measure work for computing the PMF. The conventional SMD method suffers from poor sampling efficiency and a high amount of dissipative work in regard to translocation of poly-Ns through the OmpF nanopore. Drawing inspiration from the D-MetaD approach,41 we present a joint MT and SMD method to improve the efficiency of sampling and decrease the amount of dissipative work (as a sign of divergence from the equilibrium state) in PMF calculation (Supporting Information, Supporting Method). The Jarzinsky equality (JE) and its augmented versions are

Figure 8. Joint MT and SMD-derived potential of mean force. The potential of mean force is computed for translocation of poly-N through OmpF as a function of the SMD atom’s distance from the COM of porin in angstroms. The ionic charge density of OmpF is indicated just before poly-N reaches OmpF’s opening. The densities of positive and negative ions around the OmpF’s pore are represented by the blue surface and red mesh, respectively. The gray error bars represent the SEM of four independent simulations.

The joint MT and SMD-based PMFs show a high energy barrier along the translocation path of each poly-N through OmpF. The computed PMFs indicate that rG senses a smaller barrier than dC and dG do. Such a discrepancy among the energy barriers is obvious when SMD atoms of poly-Ns reach the inner region of OmpF (error bars with no overlaps). The presence of a positive EPS around the outer region of OmpF and the attraction of a cloud of negative ions by the pore suggest that there is a probable site in the pore with affinity for poly-N. However, such a region is not detected clearly by joint MT and SMD-based PMF. Our information about the PMF of the pre-CZ region seems scarce. The adaptive biasing force (ABF) method is a versatile and powerful approach to estimating the PMF of molecules that translocate through a predetermined path such as poly-Ns that traverse a nanopore.40,92−94 Using ABF, we reasonably alleviate the sampling problem for computing an accurate PMF in a condition close to equilibrium without using the JE. In the ABF method, the system force is calculated in a predetermined direction. After specific time steps, an external biasing force is applied to neutralize the instantaneous force. The biasing force is used to compute the PMF of the process (see methods for details of the ABF). The ABF-calculated PMF of poly-Ns’ translocation through OmpF is represented in Figure 9. The total energy barriers for translocation of poly-Ns that are derived from the calculated PMFs are similar in the ABF and joint MT and SMD approach. In addition to the smaller error of ABF-based PMFs than that of joint MT and SMD ones, the obvious difference between joint MT and SMD-derived PMFs and ABF-based counterparts is the appearance of binding sites with high affinities for poly-Ns in ABF-based PMFs. Such binding sites reside in the pre-CZ region, as inferred from the EPS of porin and the distribution of ions around the protein. Poly-Ns are differentiated by two OmpF distinct regions: the outer region of OmpF that discriminates organic bases and the inner region that differentiates poly-Ns’ backbones. The K

DOI: 10.1021/acs.jpcb.5b00763 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

(Figure 10, table of ΔG_in). We remember that, because of the flexibility of the outer region’s loops, poly-Ns located at pre-CZ do not block the OmpF’s opening completely (Figure 5). For a complete translocation, poly-Ns should escape the pre-CZ trap and pass through the CZ. The poly-N molecule attracted to pre-CZ encounters a high free energy barrier. The ΔG_out metric represents the free energy barrier to poly-N passing through the CZ (Figure 10). The easiness of this step is determined by the backbone of poly-Ns; deoxy-ribose nucleic acids encounter a higher barrier than rG does. ΔG_net determines the total height of the free energy barrier to a complete translocation. The least amounts of external driving forces needed for complete translocation are estimated on the basis of ΔG_net free energies (Figure 10, table of ΔG_net). The ABF computation indicates that dC requires the highest amount of driving force, whereas rG needs the lowest amount of external driving force for translocation. In 2012, Mayer and co-workers indicated that the λ phage may find OmpF as an alternative receptor for binding to bacteria.95 We infer from their works that OmpF porin or OmpF in association with a complex may act as a way for DNA injection to periplasmic space of bacterial membrane. On the other hand, it was indicated that λ phage injects its DNA with high speed.96 These observations imply that if sufficient external forces are available DNA can pass through OmpF. Our computed PMFs for pentameric poly-Ns indicate that high amounts of driving forces are required to push DNA in OmpF’s lumen and the injection force of λ phage alone seems insufficient for translocation of DNA. We think that an auxiliary force is required. The huge energy barrier to translocation of poly-Ns appears in the CZ region and continues for about 30 Å. The rG

Figure 9. ABF-based potential of mean force. PMFs of poly-Ns’ translocation through the lumen of OmpF are computed by the ABF method. The PMF values are represented as a function of the distance of poly-N’s COM from OmpF’s COM in angstroms. Snapshots of the poly-N−protein complex during the translocation process of each type of poly-N are depicted in the inset in which poly-N is located before the outer region, at CZ, and after the inner region of OmpF from left to right. OmpF is represented by cartoon style. The bases and backbone of poly-N are represented by blue and yellow surfaces, respectively. The gray error bars represent the standard error based on eq 5 in the methods. The PMFs are smoothed by the adjacent averaging method.

thermodynamics of poly-Ns’ binding to OmpF is summarized in Figure 10. Poly-Ns are attracted to the pre-CZ region below the external loops of OmpF. The ΔG_in values indicate the free energies of binding of poly-Ns to pre-CZ. These values suggest that the binding of poly-Ns to OmpF’s entrance is favorable

Figure 10. Thermodynamics of poly-Ns’ translocation. The ΔG values of moving toward CZ (ΔG_in) and crawling from CZ toward the cytoplasm (ΔG_out) are presented beside blue arrows. These values are obtained from ABF-based PMF curves. The overall free energy barrier of the translocation process (ΔG_net) is presented on top of the orange arrow together with the minimal external driving force required for forcing polyN’s through porin’s lumen. OmpF porin is represented by the black meshed surface, and pre-CZ (right) and CZ (left) regions are highlighted by yellow solid surfaces. The error computed using eq 5 for ΔG_net is less than 30 kcal/mol, and for the other two types, ΔG is less than 20 kcal/mol. L

DOI: 10.1021/acs.jpcb.5b00763 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

and sequence of poly-N. Besides, if poly-N moves from the inner to the outer region of OmpF or if it enters OmpF from the 3′ end instead of the 5′ end, we expect a different height of the energy barrier ahead of translocation (as we are observing in a primary study). Also, different PMFs are expected for heteropolymers that may be difficult to be extrapolated from homopolymer data.

molecule encounters the smallest barrier to pass the lumen of OmpF. To supply a sufficient driving force to overcome the energy barrier ahead of rG, an estimated electrical potential of about 1 V is needed. Such voltage is higher than the threshold of stability for common biomembranes in the lab. These calculations suggest that low voltages are not able to force homopentameric poly-Ns through the OmpF channel. To test the difficulty of poly-Ns’ translocation through OmpF, we examine the possibility of translocation of dC by applying a 2 V external electric field across the simulation box. After pushing dC for 10 ns, it passes the OmpF’s CZ by tearing the OmpF structure and adjacent membrane (which was predictable). Although the molecular dynamics-based results are not similar to experiments quantitatively, it seems that MD results are not far away from the borderlines of experimental confidence intervals. dC decamers blocked the OmpF channel more than dG counterparts did in experiment.33 Our observation in which the dC poly-N needs a high applied voltage for a complete translocation may rationalize the experimental observations. dC shows the highest ΔG_net and therefore needs the greatest driving force for translocation among the studied poly-Ns (Figure 10). Consequently, dC should remain longer than other poly-Ns do at CZ when the same external voltage is applied to force poly-Ns through OmpF.

■

ASSOCIATED CONTENT

S Supporting Information *

Supporting method: The joint MT and SMD approach is describe in detail. Figure S1: The distribution of OmpF conductance in the presence of poly-Ns. Figure S2: EDA-based cross-correlation heat map of OmpF’s structure in the presence and absence of poly-Ns. Figure S3: Pair interaction energies of nucleotide MRs. Figure S4: Variations in base stacking. Figure S5: The CAVER-computed tunnels of OmpF. Figure S6: The electrostatic potential surfaces of OmpF. Figure S7: Comparison of the searching efficiency of SMD, zMT, and joint MT and SMD. Figure S8: The distributions of external work values derived from SMD, zMT, and joint MT and SMD for poly-Ns’ translocation process. Animation: The translocation of dG through OmpF porin. The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.5b00763.

■

■

CONCLUSIONS Nanopores open a new way to fast and low-cost sequencing of poly-Ns. Biomembrane channels and pores are becoming useful tools in the field of poly-N sequencing. Computational studies help us to consider the poly-N−nanopore interactions in atomic detail. There is an experimental report that suggests use of the OmpF porin in DNA sequencing.33 Here, we study the consequences of homopentameric poly-Ns’ translocation through OmpF by computational approaches for the first time. The dynamical and structural aspects of possible translocation of homopentameric poly-Ns through OmpF are studied in detail. We explore the dynamics, thermodynamics, and structural features of the translocation process to find some poly-N-specific aspects during translocation. The multidimensional approaches of the current study detect poly-N-specific behavior of OmpF during the translocation of poly-N. The PMFs of the translocation process indicate that large energy barriers exist ahead of poly-Ns. The height of computed energy barriers to poly-Ns’ translocation implies a little likelihood of translocation of homopentameric nucleic acids through OmpF by applying small external voltages. However, OmpF may still be a suitable candidate for sequencing nucleic acids. For example, if we assume that the end of a DNA molecule is attached to an AFM’s cantilever and pulled through an OmpF embedded in synthetic membrane, measuring the conductance or the required forces during pulling DNA through OmpF’s lumen may nominate OmpF for poly-N sequencing. Scientists use various methods to decrease the speed of poly-Ns’ translocation through nanopores, whereas the OmpF porin shows a considerable resistance against passing pentameric poly-Ns. In addition, the introduced MRs in the current study are options for directed engineering of OmpF to increase its poly-N specificity or boost the translocation speed. We should note that the current study depicts a special case of poly-Ns’ translocation through OmpF porin. The translocation process of poly-N is probably affected by the length

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Phone: +982433153316. Fax: +982433153322. Author Contributions

The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS Authors acknowledge the Institute for Advanced Studies in Basic Sciences (IASBS) at Zanjan, Iran, for supporting this work.

■

M

ABBREVIATIONS ABF, adaptive biasing force APC, affinity propagation clustering BFE, binding free energy COM, center of mass CZ, constriction zone dC, single-stranded poly-C DNA dG, single-stranded poly-G DNA EDA, essential dynamics analysis EPS, electrostatic potential surface JE, Jarzinsky equality LIE, linear interaction energy MR, master residue MT, metadynamics PCA, principal component analysis PME, particle-mesh Ewald PMF, potential of mean forces Poly-N, poly-nucleotide PSN, protein structure network rG, single-stranded poly-G RNA DOI: 10.1021/acs.jpcb.5b00763 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

(20) Inokuchi, K.; Mutoh, N.; Matsuyama, S.; Mizushima, S. Primary structure of the ompF gene that codes for a major outer membrane protein of Escherichia coli K-12. Nucleic Acids Res. 1982, 10, 6957− 6968. (21) Cowan, S. W.; Garavito, R. M.; Jansonius, J. N.; Jenkins, J. A.; Karlsson, R.; Konig, N.; Pai, E. F.; Pauptit, R. A.; Rizkallah, P. J.; Rosenbusch, J. P.; et al. The structure of OmpF porin in a tetragonal crystal form. Structure 1995, 3, 1041−1050. (22) Simonet, V.; Mallea, M.; Pages, J. M. Substitutions in the eyelet region disrupt cefepime diffusion through the Escherichia coli OmpF channel. Antimicrob. Agents Chemother. 2000, 44, 311−315. (23) Zhao, Z. P.; Liu, T. T.; Zhang, L.; Luo, M.; Nie, X.; Li, Z. X.; Pan, Y. High-grade mutant OmpF induces decreased bacterial survival rate. Acta Biochim. Polym. 2014, 61, 369−373. (24) Zakharov, S. D.; Eroukova, V. Y.; Rokitskaya, T. I.; Zhalnina, M. V.; Sharma, O.; Loll, P. J.; Zgurskaya, H. I.; Antonenko, Y. N.; Cramer, W. A. Colicin occlusion of OmpF and TolC channels: outer membrane translocons for colicin import. Biophys. J. 2004, 87, 3901−3911. (25) Zakharov, S. D.; Sharma, O.; Zhalnina, M.; Yamashita, E.; Cramer, W. A. Pathways of colicin import: utilization of BtuB, OmpF porin and the TolC drug-export protein. Biochem. Soc. Trans. 2012, 40, 1463−1468. (26) Kumar, A.; Hajjar, E.; Ruggerone, P.; Ceccarelli, M. Molecular simulations reveal the mechanism and the determinants for ampicillin translocation through OmpF. J. Phys. Chem. B 2010, 114, 9608−9616. (27) Mahendran, K. R.; Hajjar, E.; Mach, T.; Lovelle, M.; Kumar, A.; Sousa, I.; Spiga, E.; Weingart, H.; Gameiro, P.; Winterhalter, M.; et al. Molecular basis of enrofloxacin translocation through OmpF, an outer membrane channel of Escherichia coli–when binding does not imply translocation. J. Phys. Chem. B 2010, 114, 5170−5179. (28) Singh, P. R.; Ceccarelli, M.; Lovelle, M.; Winterhalter, M.; Mahendran, K. R. Antibiotic permeation across the OmpF channel: Modulation of the affinity site in the presence of magnesium. J. Phys. Chem. B 2012, 116, 4433−4438. (29) Robertson, K. M.; Tieleman, D. P. Orientation and interactions of dipolar molecules during transport through OmpF porin. FEBS Lett. 2002, 528, 53−57. (30) Ziervogel, B. K.; Roux, B. The binding of antibiotics in OmpF porin. Structure 2013, 21, 76−87. (31) Dhakshnamoorthy, B.; Ziervogel, B. K.; Blachowicz, L.; Roux, B. A structural study of ion permeation in OmpF porin from anomalous X-ray diffraction and molecular dynamics simulations. J. Am. Chem. Soc. 2013, 135, 16561−16568. (32) Alcaraz, A.; Nestorovich, E. M.; Lopez, M. L.; Garcia-Gimenez, E.; Bezrukov, S. M.; Aguilella, V. M. Diffusion, exclusion, and specific binding in a large channel: A study of OmpF selectivity inversion. Biophys. J. 2009, 96, 56−66. (33) Hadi-Alijanvand, S.; Mobasheri, H.; Hadi-Alijanvand, H. Application of OmpF nanochannel forming protein in polynucleotide sequence recognition. J. Mol. Recognit. 2014, 27, 575−587. (34) Jo, S.; Lim, J. B.; Klauda, J. B.; Im, W. CHARMM-GUI Membrane Builder for mixed bilayers and its application to yeast membranes. Biophys. J. 2009, 97, 50−58. (35) Wu, E. L.; Cheng, X.; Jo, S.; Rui, H.; Song, K. C.; DavilaContreras, E. M.; Qi, Y.; Lee, J.; Monje-Galvan, V.; Venable, R. M.; et al. CHARMM-GUI Membrane Builder toward realistic biological membrane simulations. J. Comput. Chem. 2014, 35, 1997−2004. (36) Jo, S.; Kim, T.; Im, W. Automated builder and database of protein/membrane complexes for molecular dynamics simulations. PLoS One 2007, 2, e880. (37) Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kale, L.; Schulten, K. Scalable molecular dynamics with NAMD. J. Comput. Chem. 2005, 26, 1781− 1802. (38) Mackerell, A. D., Jr.; Feig, M.; Brooks, C. L., 3rd. Extending the treatment of backbone energetics in protein force fields: Limitations of gas-phase quantum mechanics in reproducing protein conformational

SEM, standard error of mean SMD, steered molecular dynamics

■

REFERENCES

(1) Bayley, H.; Cremer, P. S. Stochastic sensors inspired by biology. Nature 2001, 413, 226−230. (2) Branton, D.; Deamer, D. W.; Marziali, A.; Bayley, H.; Benner, S. A.; Butler, T.; Di Ventra, M.; Garaj, S.; Hibbs, A.; Huang, X.; et al. The potential and challenges of nanopore sequencing. Nat. Biotechnol. 2008, 26, 1146−1153. (3) Hirschhorn, J. N. Genomewide association studies–illuminating biologic pathways. N. Engl. J. Med. 2009, 360, 1699−1701. (4) Venkatesan, B. M.; Bashir, R. Nanopore sensors for nucleic acid analysis. Nat. Nanotechnol. 2011, 6, 615−624. (5) Schreiber, J.; Wescoe, Z. L.; Abu-Shumays, R.; Vivian, J. T.; Baatar, B.; Karplus, K.; Akeson, M. Error rates for nanopore discrimination among cytosine, methylcytosine, and hydroxymethylcytosine along individual DNA strands. Proc. Natl. Acad. Sci. U. S. A. 2013, 110, 18910−18915. (6) Stoddart, D.; Maglia, G.; Mikhailova, E.; Heron, A. J.; Bayley, H. Multiple base-recognition sites in a biological nanopore: two heads are better than one. Angew. Chem., Int. Ed. Engl. 2010, 49, 556−559. (7) Kasianowicz, J. J.; Brandin, E.; Branton, D.; Deamer, D. W. Characterization of individual polynucleotide molecules using a membrane channel. Proc. Natl. Acad. Sci. U. S. A. 1996, 93, 13770− 13773. (8) Deamer, D. W.; Akeson, M. Nanopores and nucleic acids: prospects for ultrarapid sequencing. Trends Biotechnol. 2000, 18, 147− 151. (9) Li, J.; Gershow, M.; Stein, D.; Brandin, E.; Golovchenko, J. A. DNA molecules and configurations in a solid-state nanopore microscope. Nat. Mater. 2003, 2, 611−615. (10) Li, J.; Stein, D.; McMullan, C.; Branton, D.; Aziz, M. J.; Golovchenko, J. A. Ion-beam sculpting at nanometre length scales. Nature 2001, 412, 166−169. (11) Venta, K.; Shemer, G.; Puster, M.; Rodriguez-Manzo, J. A.; Balan, A.; Rosenstein, J. K.; Shepard, K.; Drndic, M. Differentiation of short, single-stranded DNA homopolymers in solid-state nanopores. ACS Nano 2013, 7, 4629−4636. (12) Min, S. K.; Kim, W. Y.; Cho, Y.; Kim, K. S. Fast DNA sequencing with a graphene-based nanochannel device. Nat. Nanotechnol. 2011, 6, 162−165. (13) Purnell, R. F.; Schmidt, J. J. Discrimination of single base substitutions in a DNA strand immobilized in a biological nanopore. ACS Nano 2009, 3, 2533−2538. (14) Butler, T. Z.; Pavlenok, M.; Derrington, I. M.; Niederweis, M.; Gundlach, J. H. Single-molecule DNA detection with an engineered MspA protein nanopore. Proc. Natl. Acad. Sci. U. S. A. 2008, 105, 20647−20652. (15) Bates, M.; Burns, M.; Meller, A. Dynamics of DNA molecules in a membrane channel probed by active control techniques. Biophys. J. 2003, 84, 2366−2372. (16) Stoddart, D.; Heron, A. J.; Mikhailova, E.; Maglia, G.; Bayley, H. Single-nucleotide discrimination in immobilized DNA oligonucleotides with a biological nanopore. Proc. Natl. Acad. Sci. U. S. A. 2009, 106, 7702−7707. (17) Ashkenasy, N.; Sanchez-Quesada, J.; Bayley, H.; Ghadiri, M. R. Recognizing a single base in an individual DNA strand: A step toward DNA sequencing in nanopores. Angew. Chem., Int. Ed. Engl. 2005, 44, 1401−1404. (18) de Zoysa, R. S.; Jayawardhana, D. A.; Zhao, Q.; Wang, D.; Armstrong, D. W.; Guan, X. Slowing DNA translocation through nanopores using a solution containing organic salts. J. Phys. Chem. B 2009, 113, 13332−13336. (19) Chen, R.; Kramer, C.; Schmidmayr, W.; Chen-Schmeisser, U.; Henning, U. Primary structure of major outer-membrane protein I (ompF protein, porin) of Escherichia coli B/r. Biochem. J. 1982, 203, 33−43. N

DOI: 10.1021/acs.jpcb.5b00763 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B distributions in molecular dynamics simulations. J. Comput. Chem. 2004, 25, 1400−1415. (39) Vanommeslaeghe, K.; Hatcher, E.; Acharya, C.; Kundu, S.; Zhong, S.; Shim, J.; Darian, E.; Guvench, O.; Lopes, P.; Vorobyov, I.; et al. CHARMM general force field: A force field for drug-like molecules compatible with the CHARMM all-atom additive biological force fields. J. Comput. Chem. 2010, 31, 671−690. (40) Martin, H. S. C.; Jha, S.; Coveney, P. V. Comparative analysis of nucleotide translocation through protein nanopores using steered molecular dynamics and an adaptive biasing force. J. Comput. Chem. 2014, 35, 692−702. (41) Moradi, M.; Tajkhorshid, E. Driven metadynamics: Reconstructing equilibrium free energies from driven adaptive-bias simulations. J. Phys. Chem. Lett. 2013, 4, 1882−1887. (42) Gullingsrud, J. R.; Braun, R.; Schulten, K. Reconstructing potentials of mean force through time series analysis of steered molecular dynamics simulations. J. Comput. Phys. 1999, 151, 190−211. (43) Hummer, G.; Szabo, A. Free energy reconstruction from nonequilibrium single-molecule pulling experiments. Proc. Natl. Acad. Sci. U. S. A. 2001, 98, 3658−3661. (44) Callenberg, K. M.; Choudhary, O. P.; de Forest, G. L.; Gohara, D. W.; Baker, N. A.; Grabe, M. APBSmem: A graphical interface for electrostatic calculations at the membrane. PLoS One 2010, 5, e12722. (45) Baker, N. A.; Sept, D.; Joseph, S.; Holst, M. J.; McCammon, J. A. Electrostatics of nanosystems: application to microtubules and the ribosome. Proc. Natl. Acad. Sci. U. S. A. 2001, 98, 10037−10041. (46) Aksimentiev, A.; Schulten, K. Imaging alpha-hemolysin with molecular dynamics: ionic conductance, osmotic permeability, and the electrostatic potential map. Biophys. J. 2005, 88, 3745−3761. (47) Hille, B. Ion channels of excitable membranes, 3rd ed.; Sinauer: Sunderland, MA, 2001. (48) Anishkin, A.; Sukharev, S. Water dynamics and dewetting transitions in the small mechanosensitive channel MscS. Biophys. J. 2004, 86, 2883−2895. (49) Chovancova, E.; Pavelka, A.; Benes, P.; Strnad, O.; Brezovsky, J.; Kozlikova, B.; Gora, A.; Sustr, V.; Klvana, M.; Medek, P.; et al. CAVER 3.0: A tool for the analysis of transport pathways in dynamic protein structures. PLoS Comput. Biol. 2012, 8, e1002708. (50) Kozlikova, B.; Sebestova, E.; Sustr, V.; Brezovsky, J.; Strnad, O.; Daniel, L.; Bednar, D.; Pavelka, A.; Manak, M.; Bezdeka, M.; et al. CAVER Analyst 1.0: Graphic tool for interactive visualization and analysis of tunnels and channels in protein structures. Bioinformatics 2014, 30, 2684−2685. (51) Comer, J.; Gumbart, J. C.; Henin, J.; Lelievre, T.; Pohorille, A.; Chipot, C. The adaptive biasing force method: everything you always wanted to know but were afraid to ask. J. Phys. Chem. B 2015, 119, 1129−1151. (52) Rodriguez-Gomez, D.; Darve, E.; Pohorille, A. Assessing the efficiency of free energy calculation methods. J. Chem. Phys. 2004, 120, 3563−3578. (53) Straatsma, T. P.; Berendsen, H. J. C.; Stam, A. J. Estimation of statistical errors in molecular simulation calculations. Mol. Phys. 1986, 57, 89−95. (54) Smith, E. B.; Wells, B. H. Estimating errors in molecular simulation calculations. Mol. Phys. 1984, 52, 701−704. (55) Almlöf, M.; Carlsson, J.; Åqvist, J. Improving the accuracy of the linear interaction energy method for solvation free energies. J. Chem. Theory Comput. 2007, 3, 2162−2175. (56) Hansson, T.; Marelius, J.; Aqvist, J. Ligand binding affinity prediction by linear interaction energy methods. J. Comput.-Aided Mol. Des. 1998, 12, 27−35. (57) Aqvist, J.; Medina, C.; Samuelsson, J. E. A new method for predicting binding affinity in computer-aided drug design. Protein Eng. 1994, 7, 385−391. (58) Bakan, A.; Dutta, A.; Mao, W.; Liu, Y.; Chennubhotla, C.; Lezon, T. R.; Bahar, I. Evol and ProDy for bridging protein sequence evolution and structural dynamics. Bioinformatics 2014, 30, 2681− 2683.

(59) Bakan, A.; Meireles, L. M.; Bahar, I. ProDy: protein dynamics inferred from theory and experiments. Bioinformatics 2011, 27, 1575− 1577. (60) Amadei, A.; Linssen, A. B.; Berendsen, H. J. Essential dynamics of proteins. Proteins 1993, 17, 412−425. (61) Hayward, S.; de Groot, B. L. Normal modes and essential dynamics. Methods Mol. Biol. 2008, 443, 89−106. (62) Lu, X. J.; Olson, W. K. 3DNA: A versatile, integrated software system for the analysis, rebuilding and visualization of threedimensional nucleic-acid structures. Nat. Protoc. 2008, 3, 1213−27. (63) Lu, X. J.; Olson, W. K. 3DNA: A software package for the analysis, rebuilding and visualization of three-dimensional nucleic acid structures. Nucleic Acids Res. 2003, 31, 5108−5121. (64) Brinda, K. V.; Vishveshwara, S. A network representation of protein structures: Implications for protein stability. Biophys. J. 2005, 89, 4159−4170. (65) Seeber, M.; Cecchini, M.; Rao, F.; Settanni, G.; Caflisch, A. Wordom: A program for efficient analysis of molecular dynamics simulations. Bioinformatics 2007, 23, 2625−2627. (66) Seeber, M.; Felline, A.; Raimondi, F.; Muff, S.; Friedman, R.; Rao, F.; Caflisch, A.; Fanelli, F. Wordom: A user-friendly program for the analysis of molecular structures, trajectories, and free energy surfaces. J. Comput. Chem. 2011, 32, 1183−1194. (67) Frey, B. J.; Dueck, D. Clustering by passing messages between data points. Science 2007, 315, 972−976. (68) Neidle, S. Principles of nucleic acid structure; Elsevier/Academic Press: Amsterdam, London, 2008. (69) Bodenhofer, U.; Kothmeier, A.; Hochreiter, S. APCluster: An R package for affinity propagation clustering. Bioinformatics 2011, 27, 2463−2464. (70) Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual molecular dynamics. J. Mol. Graphics 1996, 14, 33−38. (71) Aksimentiev, A. Deciphering ionic current signatures of DNA transport through a nanopore. Nanoscale 2010, 2, 468−483. (72) Lu, H.; Schulten, K. Steered molecular dynamics simulations of force-induced protein domain unfolding. Proteins 1999, 35, 453−463. (73) Ribeiro, R. F.; Marenich, A. V.; Cramer, C. J.; Truhlar, D. G. The solvation, partitioning, hydrogen bonding, and dimerization of nucleotide bases: A multifaceted challenge for quantum chemistry. Phys. Chem. Chem. Phys. 2011, 13, 10908−10922. (74) Muthukumar, M. Theory of capture rate in polymer translocation. J. Chem. Phys. 2010, 132, 195101. (75) Bonthuis, D. J.; Zhang, J.; Hornblower, B.; Mathe, J.; Shklovskii, B. I.; Meller, A. Self-energy-limited ion transport in subnanometer channels. Phys. Rev. Lett. 2006, 97, 128104. (76) Mathe, J.; Aksimentiev, A.; Nelson, D. R.; Schulten, K.; Meller, A. Orientation discrimination of single-stranded DNA inside the alphahemolysin membrane channel. Proc. Natl. Acad. Sci. U. S. A. 2005, 102, 12377−12382. (77) Basle, A.; Iyer, R.; Delcour, A. H. Subconductance states in OmpF gating. Biochim. Biophys. Acta 2004, 1664, 100−107. (78) Benz, R. Bacterial and eukaryotic porins: Structure, function, mechanism; Wiley-VCH: Weinheim, Germany, 2004. (79) Tomita, N.; Mohammad, M. M.; Niedzwiecki, D. J.; Ohta, M.; Movileanu, L. Does the lipid environment impact the open-state conductance of an engineered beta-barrel protein nanopore? Biochim. Biophys. Acta 2013, 1828, 1057−1065. (80) Alcaraz, A.; Queralt-Martin, M.; Garcia-Gimenez, E.; Aguilella, V. M. Increased salt concentration promotes competitive block of OmpF channel by protons. Biochim. Biophys. Acta 2012, 1818, 2777− 2782. (81) Markosyan, S.; De Biase, P. M.; Czapla, L.; Samoylova, O.; Singh, G.; Cuervo, J.; Tieleman, D. P.; Noskov, S. Y. Effect of confinement on DNA, solvent and counterion dynamics in a model biological nanopore. Nanoscale 2014, 6, 9006−9016. (82) Aguilella-Arzo, M.; Andrio, A.; Aguilella, V. M.; Alcaraz, A. Dielectric saturation of water in a membrane protein channel. Phys. Chem. Chem. Phys. 2009, 11, 358−365. O

DOI: 10.1021/acs.jpcb.5b00763 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B (83) Lindorff-Larsen, K.; Piana, S.; Dror, R. O.; Shaw, D. E. How fast-folding proteins fold. Science 2011, 334, 517−520. (84) Hadi-Alijanvand, H.; Proctor, E. A.; Goliaei, B.; Dokholyan, N. V.; Moosavi-Movahedi, A. A. Thermal unfolding pathway of PHD2 catalytic domain in three different PHD2 species: computational approaches. PLoS One 2012, 7, e47061. (85) Shanahan, H. P.; Garcia, M. A.; Jones, S.; Thornton, J. M. Identifying DNA-binding proteins using structural motifs and the electrostatic potential. Nucleic Acids Res. 2004, 32, 4732−4741. (86) Guenot, J.; Fletterick, R. J.; Kollman, P. A. A negative electrostatic determinant mediates the association between the Escherichia coli trp repressor and its operator DNA. Protein Sci. 1994, 3, 1276−1285. (87) Harris, R. C.; Mackoy, T.; Dantas Machado, A. C.; Xu, D.; Rohs, R.; Fenley, M. O. Chapter 3 Opposites attract: Shape and electrostatic complementarity in protein-DNA complexes. Innovations in Biomolecular Modeling and Simulations: Vol. 2; The Royal Society of Chemistry: Cambridge, U.K., 2012; Vol. 2, pp 53−80. (88) Minh, D. D.; McCammon, J. A. Springs and speeds in free energy reconstruction from irreversible single-molecule pulling experiments. J. Phys. Chem. B 2008, 112, 5892−5897. (89) Park, S.; Schulten, K. Calculating potentials of mean force from steered molecular dynamics simulations. J. Chem. Phys. 2004, 120, 5946−5961. (90) Kosztin, I.; Barz, B.; Janosi, L. Calculating potentials of mean force and diffusion coefficients from nonequilibrium processes without Jarzynski’s equality. J. Chem. Phys. 2006, 124, 064106. (91) Holland, B. W.; Gray, C. G.; Tomberli, B. Calculating diffusion and permeability coefficients with the oscillating forward-reverse method. Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys. 2012, 86, 036707. (92) Darve, E.; Rodriguez-Gomez, D.; Pohorille, A. Adaptive biasing force method for scalar and vector free energy calculations. J. Chem. Phys. 2008, 128, 144120. (93) Hénin, J.; Fiorin, G.; Chipot, C.; Klein, M. L. Exploring multidimensional free energy landscapes using time-dependent biases on collective variables. J. Chem. Theory Comput. 2009, 6, 35−47. (94) Dehez, F.; Pebay-Peyroula, E.; Chipot, C. Binding of ADP in the mitochondrial ADP/ATP carrier is driven by an electrostatic funnel. J. Am. Chem. Soc. 2008, 130, 12725−12733. (95) Meyer, J. R.; Dobias, D. T.; Weitz, J. S.; Barrick, J. E.; Quick, R. T.; Lenski, R. E. Repeatability and contingency in the evolution of a key innovation in phage lambda. Science 2012, 335, 428−432. (96) Grayson, P.; Han, L.; Winther, T.; Phillips, R. Real-time observations of single bacteriophage lambda DNA ejections in vitro. Proc. Natl. Acad. Sci. U. S. A. 2007, 104, 14652−14657.

P

DOI: 10.1021/acs.jpcb.5b00763 J. Phys. Chem. B XXXX, XXX, XXX−XXX