Weak Selectivity Predicted for Modeled Bundles of ... - ACS Publications

Dec 1, 2016 - Dhani Ram Mahato and Wolfgang B. Fischer*. Institute of Biophotonics and Biophotonics & Molecular Imaging Research Center (BMIRC), ...
1 downloads 0 Views 3MB Size
Subscriber access provided by Olson Library | Northern Michigan University

Article

Weak Selectivity Predicted for Modeled Bundles of the Viral Channel Forming Protein E5 of Human Papillomavirus-16 Dhani Ram Mahato, and Wolfgang B. Fischer J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.6b10050 • Publication Date (Web): 01 Dec 2016 Downloaded from http://pubs.acs.org on December 2, 2016

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry B is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Weak Selectivity Predicted for Modeled Bundles of the Viral Channel Forming Protein E5 of Human Papillomavirus-16.

Dhani Ram Mahato and Wolfgang B. Fischer*

Institute of Biophotonics and Biophotonics & Molecular Imaging Research Center (BMIRC), School of Biomedical Science and Engineering, National Yang-Ming University, Taipei, Taiwan

*Correspondence to: W. B. Fischer, Institute of Biophotonics, School of Biomedical Science and Engineering, National Yang-Ming University, 155, Li-Nong St., Sec. 2, Taipei, 112, Taiwan

E-mail address: [email protected] Tel: +886-2-2826-7394, Fax: +886-2-2823-5460

ACS Paragon Plus Environment

1

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 38

Abstract: Protein E5 is a polytopic 83 amino acid membrane protein with three transmembrane domains (TMDs) encoded by the high-risk human papillomavirus 16 (HPV-16). HPV-16 is found to be the causative agent for cervical cancer. Protein E5 is, amongst other proteins (e.g. E6, E7), expressed at an ‘early (E) stage when the cell turns malignant. It has been experimentally found that E5 forms hexameric assemblies which show the characteristics of the class of so called channel forming proteins by rendering lipid membranes permeable to ions and small molecules. Protein E5 is used to achieve structural models of the protein in assembled bundles using a force field based docking approach. Extended MD simulations of selected bundles in fully hydrated lipid bilayers suggest the second TMD to be pore lining allowing for water columns to exist within the lumen of the pore. Full correlation analysis (FCA) suggests asymmetric dynamics within the monomers of the bundle. Potential of mean force (PMF) calculations of a snapshot structure of the putative open pore of the protein bundle suggest low selectivity.

ACS Paragon Plus Environment

2

Page 3 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Introduction

Ion channels are integral membrane proteins which enable the passive flux of ions across the lipid membranes (see

1

and references therein). A series of computational studies report about the

contributions of biophysical parameters, such as pore size, electrostatics or hydrophobicity within the pore, to ion flux

2-5

. In addition, molecular dynamics (MD) simulations contribute to the

understanding of the mechanism of function by exploring the dynamics of the proteins involved in ion translocation and selectivity 6-7.

HPV-16 E5 is an 83 amino acid oncoprotein encoded by human papilloma virus type 16 (HPV-16) 8 and belongs to the class of viral channel forming proteins

6-7, 9

HPV-16 is the main cause of

cervical cancer (for a review and treatment see 10-11) and in an early stage (abbreviated as E) of viral infection three proteins are expressed, E5, E6 and E7. E5 is found in the endoplasmic reticulum and also suggested to associate with the Golgi apparatus

12

and the perinuclear region

13

Three

hydrophobic regions are identified in E5 from hydropathy index calculations of the amino acids sequence of E5 8, 14, one longer stretch of residues 8-30, and two stretches of residues 37-52 and 5876, eventually too short to stretch the membrane

15

. The third hydrophobic region is proposed to

contain a voltage gating motif similar to the region M2 in connexins 16-17.

Its orientation is determined to present the C terminus to the cytoplasm and the N terminus to the lumen of the ER

13

. From native PAGE analysis E5 is found to adopt dimeric and hexameric

assembly with the latter forming defined circular arrangement visible in transmission electron microscopy (TEM) images 9.

ACS Paragon Plus Environment

3

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 38

The structural features of the TMDs of E5 have reported to be helical based on FTIR spectroscopic investigations 16. This finding has been flanked by the fact that Kyte-Doolittle hydrophobicity plots reveal three hydrophobic stretches in the sequence of E5 18. In the contrary, from CD-spectroscopic investigations, a large β-sheet content for the first and third TMD is reported 19. Based on secondary structure prediction programs and the circular arrangement of the TEM data an in silico model of E5 as a hexamer in which the second TMD lining the putative pore has been reported recently 9. That model is used to predict binding sites of adamantanes and another active inhibitor all of which are identified to reduce the efflux of fluorescence dyes from E5 containing liposomes.

Expressed in an early stage of cell transformation/infection the viral DNA of E5 is often not integrated into the host genome later 20. This is seen as that E5 seems not to be necessarily essential for the malignant progression of the infected cell 15. In the early stage of malignant transformation, E5 mRNA is readily identified in HPV lesions

21

. Consequently, E5 is a potential drug target for

preventive therapy.

In this in silico study, a complete computational structural build-up

22

of E5 is presented followed

by an analysis of the characteristics of the protein in respect of its internal dynamics using full correlation analysis (FCA)

23

. Ion selectivity of an open form of the channel forming bundle is

investigated using potential of mean force (PMF) calculations and MD simulations under applied physiologically relevant voltages.

ACS Paragon Plus Environment

4

Page 5 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Materials and Methods:

The amino acid sequence of protein E5 has been taken according to the sequence in the NCBI database, Accession Number AAD33257.1 24-26:

MTNLDTASTT10

LLACFLLCFC20

VLLCVCLLIR30

PLLLSVSTYT40

SLILLVLLLW50

ITAASAFRCF60 IVYIVFVYIP70 LFLIHTHARF80 LIT

The TMDs were identified using secondary structure prediction programs (SSPPs): TMHMM DAS

28

, PSIPRED

29

, MINNOU

30

, PHOBIUS

31

27

,

, TOPO2 (Johns, S. J. (2005), TOPO2,

Transmembrane proteins display software. http://www.sacs.ucsf.edu/TOPO2), OCTOPUS

32

,

TMMOD 33 and HMMTop 34. Regions of the three TMDs were chosen as follows:

TMD18-29: STT10 LLACFLLCFC20 VLLCVCLLI, TMD235-53: SVSTYT40 SLIILVLLLW50 ITA, and TMD358-80: CF60 IVYIIFVYIP70 LFLIHTHARF80 LIT.

Software BETApro 35 and BeST 36 have been used to predict β-sheet motifs.

Helices of each of the three TMDs were generated with backbone dihedrals of ϕ = −65° and φ = −39° (ideal helix) using the integrated protein builder program of the MOE Suite (v2012.10, Chemical Computing Group, Japan). The respective residues connecting the TMDs (loop1: residues

ACS Paragon Plus Environment

5

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 38

30-34; loop2: residues 54-58) were generated onto the helices of the monomers using the program LOOPY 37.

Assembly of monomers and hexamers

The three TMDs were assembled (i) in a concerted way

38

, in which they were placed around a

central axis and the conformational space screened simultaneously with three degrees of freedom, distance, rotational angle and tilt, or (ii) in a sequential way

39-40

. In the latter case two TMDs are

assembled first followed by the assembly of the third TMD onto the dimeric TMD. In all the assembly protocols, distance between the helices was screened in steps of 0.5 Å (monomer) and 0.25 Å (hexamer) covering 7.5 Å to 14.0 Å for monomers and 12.0 Å to 20.0 Å for hexameric bundles with a step width of 0.5 Å. Rotational angles around the helical axis were altered in steps of 5° covering 360° and the tilt was altered with a step width of 4° ranging −36 ° to +36 °. Each new position was achieved by moving the Cα atoms of the helices first and reshaping the side-chains into the structures using most-likely conformations of the side-chains from an internal library of the MOE software (www.chemcomp.com). A short minimization of 15 steps of steepest decent followed in order to remove overlap of close contacts. The potential energy of each of the conformer was calculated using AMBER94 force field (ff) with a dielectric of ε = 2. Screening the conformational space of the hexameric bundles was done by moving the six monomers along the three degrees of freedom around a central six-fold pseudo symmetry axis in a concerted way.

In the protocols including MD simulations, the structure of the last frame of the respective simulation was taken forward to the next step of the protocol.

ACS Paragon Plus Environment

6

Page 7 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

N- and C-terminal segments were modeled using SABLE 41 and Jpred3 servers

42

and added to the

bundles while the hexamer was formed using Protein Builder program of MOE.

MD simulations

The monomers and bundles were embedded into fully hydrated patch of 16:0−18:1 diester PC, 1palmitory-2-oleoyl-sn-glycero-3-phospho-chloine (POPC) lipids removing lipids which overlapped with the peptides. Similar numbers of lipid molecules were taken out from both leaflets of the lipid bilayer. Lipid parameters were adopted from Chandrasekhar et al. 43.

After insertion of the protein structures into the lipid membrane (MOE suit) and removing bad overlaps if necessary, an equilibration protocol with stepwise energy minimization of 0.2 fs was adapted. Before minimization Cl-ions were added as necessary to neutralize the simulation box. The system was heated gradually from 0 K, to 250, and 310 K and kept there for 500 ps each. Then in subsequent four steps of equilibration each for the duration of 500 ns, restraining is released by applying a harmonic force of magnitude 1000, 500, 250, and 0 kJ/mol/nm, sequentially, along each of the three axis. The individual steps of this equilibration protocol sum up to 3.5 ns. All simulations, including the 100 ns production run were conducted with GROMACS 4.6.1 (University of Groningen, Netherland) using the Gromos96 ffG45a3 force field to obtain a relaxation of the conformation of the protein and the environment as much as possible. The system (NPT ensemble) was moved at a step size of 0.2 fs. Surface tension of 750 bar⋅nm was applied in xy plane whereas pressure of 1 bar was applied in z-direction. Compressibility of 4.5 × 10-5 bar-1 was

ACS Paragon Plus Environment

7

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 38

applied in all directions. A Berendsen thermostat with a coupling constant of 0.1 ps and a coupling time of 0.2 fs was used as well as the Nosé-Hoover barostat.

The monomer/hexamer systems used in the simulations consisted of 749/4494 bundle atoms, 272/412 POPC, and 8745/14602 simple point charge (SPC) water molecules, including 3/18 Clions. The total number of atoms add up to 41131/69742 atoms in the simulation boxes.

Calculation of PMF

Na-, K-, Ca-, and Cl-ions were individually pulled through the central pore of the bundles along the z-axis by using steered MD. Ions were pulled for 1000 ps with pulling force of 1000 kJ/mol⋅nm-2 and pulling rate of 5 nm/ns (50 Å/ns) to generate respective starting configurations. Frames (~100 frames) separated by steps of 0.025 nm were chosen for the umbrella sampling (US) simulations. Each of the US simulations was run for 2 ns. The collected data of these simulations were used to calculate the PMF using g_wham algorithm of the GROMACS suit. Further, bootstrap analysis (n=200) was used to estimate the statistical error.

Application of Voltage

Additional forces were applied along the z-axis to a hexamer system, hydrated by a 1 molar solution of NaCl. Applied forces corresponded to 5, 11, 33, 55, 77, 99 mV and were applied for a production run of 200 ns. g_count and g_flux programs of the GROMACS suit were used to calculate the ion permeation in the system.

ACS Paragon Plus Environment

8

Page 9 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Full correlation analysis

Full Correlation Analysis (FCA) was performed to identify functional correlated motions of the protein using the GROMACs suit as reported earlier 44. In brief, the g_fca program which uses PCA eigenvectors for the initial guess, was used for the FCA analysis of the backbone atoms

23

. The

g_covar tool was used to calculate the covariance matrix and by diagonalization of the matrix, eigenvectors and eigenvalues are obtained and ranked by size. These data are later taken as input for the initial guess for orthonormal coordinate transformation. The tool g_anaeig is applied to sort the essential frames of the trajectory along any of the eigenvectors. By performing FCA on the first 100 PCA eigenvectors, these data were ranked by the fluctuation amplitude. A minimized mutual information (MI) matrix of the order (NxN) is calculated. In this matrix N corresponds to the number of FCA vectors which were similar to the number of minimized principle components.

Hard and software

The equilibration runs were generated on a Dell Optiplex workstation with 8 core. The consequent production runs were sent to an ALPS - Acer AR585 F1 cluster (National Center for Highperformance Computing (NCHC)). Radii of the pores of the bundles were calculated using the program HOLE

45

. Plots and pictures were generated using Xmgrace (Weizmann Institute of

Science), VMD-1.9.2 (University of Illinois) and MOE 2014.09. The Figures were generated with the PyMOL program (DeLano Scientific LLC, San Francisco, CA).

ACS Paragon Plus Environment

9

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 38

Results

Structural features

The primary structure of E5 is proposed to harbor three helical segments by a series of secondary structure prediction programs (Figure 1A). Six out of nine programs assign a gap between helical regions at positions 30 – 35. Seven programs assign a second gap at position 55 – 60. Both gaps are proposed to exist of more than 2 residues. Based on the data TMDs are suggested and used in this study ranging from residues 8 – 29, 35 – 53 and 59 – 80. A β-sheet structure is not proposed by programs designed to predict this particular structural element .

Two ideal helix assemblies have been generated first, before the third helix is added (see protocols P1 and P4 as well as P2 and P5 in Table 1). In addition, a ‘concerted’ assembly in which all three TMDs approach each other simultaneously is also considered (P3 and P6). Alternatively, each of the TMDs is simulated in a 100 ns MD simulation prior to conducting the same docking protocols mentioned above (P1 – P3). Protocols P4 – P6 include ideal helices to generate the monomers. Every monomer of each assembly protocol is then copied six times and screened for the best configuration according to the ‘concerted’ docking protocol. In some of the protocols, marked as ‘a’ in Table 1, the monomers have undergone a 100 ns MD simulation before applying the docking protocol. In protocols marked as ‘b’ the monomers are used without being applied in an MD simulation.

ACS Paragon Plus Environment

10

Page 11 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

With the docking protocols the conformation space is screened in 2D applying three degrees of freedom, distance, tilt and rotational angle. The potential energy of each conformation is calculated (see Figure S1A as an example for applying P1a and P1b, as well as P4b). Due to the fact that it is relatively hydrophilic the focus is on achieving bundles in which TMD2 is pore lining which is achieved with e.g. P1a (Figure 1B) and P1b, as well as P4b. Bundles with TMD2 lining the pore are still amongst the best 12 structures generated by all the protocols (Table 1).

Bundles generated using P1a (bundle-P1a) and P1b (bundle-P1b) as well as the respective ideal bundle generated using P4b (bundle-P4b) are chosen for further structural and functional evaluation. Calculation of the RMSD values of bundles-1a/b from a 100 ns MD simulation show that values for both simulations reach a plateau at about 50 ns (Figure S1B). The formation of a compact bundle-4b structure has possibly just finished.

ACS Paragon Plus Environment

11

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 38

Fig. 1: (A) Identification of helical motifs in the sequence of E5 protein using a series of secondary structure prediction programs for α-helical motifs with a TMD, TMHMM, DASTM, PSIPRED, MINNOU, PHOBIUS, TOPO2, OCTOPUS, TMMOD and HMMTOP. Predicted residues to adopt a helical motif are marked with asterisk and an orange bar. Programs BETApro and BeST are used for predicting β-fold. (B) Structural model according to P1a (see Table 1) with TMD2 facing the pore. Ideal helices TMD1 and TMD2 are docked first. TMD3 is then docked onto the first two TMDs to generate the monomer. The monomer is then copied six times and the bundle formed by a concerted docking approach. Hydrophilic residues Ser-35, -37, -41 and -52 (red) as well as Thr-38 and -49 (yellow) are highlighted. All other residues of TMD2 are highlighted in orange. Atoms are depicted in van der Waals representation and two monomers are removed to view inside the lumen of the pore.

ACS Paragon Plus Environment

12

Page 13 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Tab. 1: Values of the lowest energy structures (Eint (internal) in kJ/mol) for each configuration derived from specific assembly protocols (P1 – P6). TMDs in squared brackets are docked first before the respective third helix is docked onto the dimeric structures. The resulting monomer (Mon) is copied 6 times to form the hexamer (Hex) and aligned separately so that each of the TMDs would align the inside of the pore (‘TMD inside’). P1 – P3 describe the docking of ideal helices into monomers, while P4 – P6 describe docking of helices which have been simulated for 100 ns. Abbreviations ‚a‘ and ‚b‘ denote that the protocol included simulations and no simulations on the monomer structure, respectively. Shaded numbers highlight the best structures, rank 1.

Mono

Hex

TMD inside

1

2

3

Eint

Rank

Eint

Rank

Eint

Rank

P1

[H1+H2]+H3

a b

-7724 -8296

2 3

-7871 -8391

1 1

-8194

9

P2

H1+[H2+H3]

a b

-7865 -7344

2 1

-7854 -

4 -

-7871 -

1 -

P3

H1+H2+H3

a b

-8779 -

6 -

-8791 -

2 -

-8887 -7872

1 1

P4

[H1+H2]+H3

a b

-9707 -6487

1 >1000

-9612 -8053

2 5

-9370 -8128

25 1

P5

H1+[H2+H3]

a b

-8327 -7019

3 17

-8131 -7825

12 2

-8446 -7843

1 1

P6

H1+H2+H3

a b

-6919 -7837

3 1

-7218 -7714

2 4

-7219 -7800

1 2

ACS Paragon Plus Environment

13

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 38

At the beginning of the 100 ns MD simulation of the three bundles, bundles-P1a and –P1b do not show any water molecules within the putative pore (Figure 2). During the 100 ns MD simulation only bundle-P1a develops a continuous water column after about 3 ns due to its large size (see Figure S2A).

Bundle-4b starts with a continuous water column at 0 ns but loses the water

molecules after about 20 ns due to the smaller pore size (Figure S2A).

Fig. 2: Hexameric bundles derived from three protocols. Bundle-P1a, bundle-P1b and bundle-P4b, are shown in a side view at the beginning (0 ns) and the end of the MD simulation (100 ns). Protein backbones are shown in yellow ribbon representation and hydrophilic residues such as serines and threonines highlighted in red stick modus. The van der Waals surface of the entire protein is shown in light grey. Water molecules are shown in blue van der Waals surface representation.

ACS Paragon Plus Environment

14

Page 15 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

FCA analysis

Visible inspection of structures belonging to the first eigenvector of a full correlation analysis (FCA) of the bundle structures done over the 100 ns MD simulations as well as over the last 20 ns of that simulation reveals that the monomers and the respective TMDs do not show collective and unidirectional movement (Figure 3 and Figure S3). Classifying the amplitude of the structures in the individual FCAs with a B-factor reveals that the averaged B-factor over all Cα-atoms and structures of eigenvector (FCA) 1 declines in the sequence: bundle-P1a > bundle-P1b > bundle-P4b. Bundle-P1a is the most flexible bundle. The same trend is observed for the data calculated for the last 20 ns. Plotting FCA2 over FCA1 for either the 100 ns trajectory or the last 20 ns trajectory shows that the structures of bundle-P1a and -P4b move in multiple smaller overlapping clusters (Figure 3). On the contrary, connecting frames of bundle-P1b reveals that localized movements of the bundle are visited consecutively (see also the values of the FCAs over time in Figure S4). Correlating the data of the B-factors with the dynamics of the protein bundle, suggests that the ‘open’ bundle-P1a is relatively dynamic while the ‘closed’ bundle is relatively rigid.

ACS Paragon Plus Environment

15

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 38

Fig. 3: Visualization of the first eigenvector of the structures calculated from a full correlation analysis (FCA), FCA1, of bundle-P1b (top row) over the entire duration of the simulations (100 ns) and the last 20 ns (80 - 100 ns) (first and last image, respectively). The color of the bundles follows the coding for B-factor calculations with blue represents low B-factor and red represents a high value. Averaged B-factors over all Cα-atoms and (i) averaged over 100 ns are: bundles-P1a ((37.3 ± 11.7) Å2), bundle-P1b ((31.8 ± 5.8) Å2), bundle-P4b ((25.2 ± 7.9) Å2), as well as (ii) averaged over the last 20 ns: bundle-P1a (18.7 ± 4.2 Å2), bundle-P1b (17.1 ± 5.8 Å2), bundleP4b (14.7 ± 2.0 Å2). Projection of the trajectory of the second FCA mode on the first FCA mode of bundle-P1a (MI = 0.35) from the entire duration of the simulation as well as the last 20 ns (MI = 0.05) (first and second plot, respectively of top row). Similar images and plots for bundle-P1b (MI = 0.44 (100 ns), MI = 0.01 (last 20 ns) (middle row) as well as bundle-P4b (MI = 0.54 (100 ns), MI = 0.00 (last 20 ns) bottom row).

ACS Paragon Plus Environment

16

Page 17 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

PMF calculations

The PMF values for ions through bundle-P1a (minimum pore radius ((0.64 ± 0.2) nm) along a linear reaction coordinate following the pore axis of the bundles show energy barriers and minima in the range of + and -0.5 kcal/mol independent of the initial pulling direction of the ions through bundle 1a (Figure 4). The PMF values of all ions overlap especially for the profile if initially pulled from the N→C terminus. All PMF curves in common is a gradual increase of the values towards the center of the bundle followed by a decrease of the values, even some of them turning negative. Correlation of the data with the structure of the pore, reveals that the hydrophobic stretch at the C terminal side of TMD2 imposes an energy barrier while the hydrophilic N terminal side imposes no barrier but rather an energetic trap for the ions.

Performing the same calculations for bundle-P4b (minimum pore radius (0.41 ± 0.13) nm) let the PMF values cover an energy range from up to +2.8 kcal/mol (e.g. Ca- and Cl-ions) to -0.75 kcal/mol (Figure 4). The overall pattern as mentioned for bundle-1a remains, an energy minimum at the C terminal side for ions being pulled from N→C and an energy barrier for the ions along the hydrophobic stretch. The monovalent cations experience a slightly larger energy barrier (~1.0 kcal/mol) when being pulled in C→N direction than for being pulled N→C (~0.75 kcal/mol). Clion is most sensitive to the hydrophilic/hydrophobic stretches when being pulled from C→N reaching maxima of ~2.75 kcal/mol. The ions experience a similar increase in energy to cross from C→N as compared for N→C, suggesting rectification behavior of the channel.

ACS Paragon Plus Environment

17

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 38

As a result, the ions seem to face lower energy barriers when approaching from the N terminal side of TMD2 and the protein bundle rather when approaching from the C terminal side. The monovalent cations show almost identical PMF profiles and are the least sensitive ions to the size of the pore. Ca- and Cl-ions are most sensitive facing larger energy barriers from the hydrophobic stretch of the N terminal side of TMD2. Pulling direction affects the PMF profile for each ion imposing an eventually rectifying conductance behavior of the bundle. In a dynamic bundle, when narrowing of the pore size this leads to an early ebbing of Ca- and Cl-ion conductance.

Fig. 4: Potential of Mean Force (PMF) calculation for Ca- (blue), Na- (red), K- (green), and Cl-ions (black) along the pore axis of bundles-P1a and -P4b. The frames for PMF calculation were taken from 20ns simulation while both the hexamers having water molecules within the entire length of the pore. The backbone of the bundles is shown in ribbon mode (white) and hydrophilic residues of TMD2 (yellow) lining the pore of hexamer are highlighted in stick mode (red). The PMF values are evaluated using bootstrap analysis (No. of bins = 200).

ACS Paragon Plus Environment

18

Page 19 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Ion dynamics

With increasing holding potentials an increasing number of Na- and Cl-ions passes through the pore, independent of the direction of the potential applied (Figure 5). Beyond an applied voltage of 55 mV the pore radius is increased but this does not lead to a further increase of the number of ions passing the pore. More Na-ions pass through the pore at negative voltage (> 40 Na-ions for > 33 mV) which corresponds to a movement from N→C than when applying positive voltage (~40 Naions for > 33 mV), corresponding to a movement from C→N (Figure 5). The number of Cl-ions is not affected by the direction of applied voltages. Application of more than 200 mV are performed leading to a continuous enlargement of the pore of the bundle during the simulation and therefor are not considered for analysis. Application of negative voltage with the Na-ions moving from N→C, the bundle develops a localized minimum (needle eye) in the middle of the pore (Figure S2B) with radii of around (0.56 ± 0.28) nm (-5 mV) to (0.11 ± 0.18) nm (-99 mV) (Table S1) caused by pore lining residues Leu-45 and Leu-49 . The mode of the transition of the ions is mostly single filed, the ions pass one after the other and only one ion is in the pore. Occasionally the two ions of opposing sign pass each other within the lumen of the pore (data not shown).

The result of the simulations is that the analysis of the bundle structure in an open state supports the notion of that hexameric E5 with TMD2 facing the pore produces a weakly selective channel. Narrowing the pore increases selectivity for monovalent cations. Eventually, the architecture of the bundle, develops a slightly rectifying channel.

ACS Paragon Plus Environment

19

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 38

Fig. 5: Number of ions permeating tough bundle-P1a during the MD simulation under applied voltages such as +/-5 mV, +/- 11 mV, +/-33 mV, +/- 55 mV, +/- 77 mV and +/- 99 mV. Bundle-P1a is embedded in a 1 M solution of NaCl (109 Na-ions / 127 Cl-ions). Positive voltage leads to the movement of the Na-ion from N→C, negative voltage leads to the movement of the Cl-ions from N→C. The view into the pore of the bundles at voltages +/- 5 mV, +/- 55 mV and +/- 99 mV are shown underneath the curve. The structures are the last frame of the 200 ns MD simulation.

ACS Paragon Plus Environment

20

Page 21 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Discussion

Assessment of the approach The structural features of the E5 protein as a monomer and its characteristics in an arrangement as a channel protein are based on a 2D docking approach based on screening of conformational space and scoring each of the structure by a force field based approach. The results are bundle models in which TMD2 is lining the pore as suggested by others9. With the generated structures ion selectivity is proposed ahead of experimental results. The docking protocol has been developed earlier

39

and

referenced against the experimentally derived structures of the TMDs of M2 of influenza A in either a closed state (PDB entry ID: 1NYJ 46) or an open state with bound channel blocker (PDB entry ID: 2H95

47

). Also when using experimentally derived structures of TMDs (TMD M2 of GLIC (PDB

entry ID: 3EHZ 48) aiming to retrieve the experimental structure has been achieved 49. Proposing an assembled oligomer of the TMD of Vpu of HIV-1 also delivers structures which can find their match with NMR spectroscopic data 50.

Based on phylogenetic analysis protein E5 can be grouped into several genera

51-52

. Protein E5 of

species HPV-16, with its e.g. genotypes 16, 31, 33, 35, 52, 58, 67, belongs to the genus alphapapillomavirus, named as E5α which comprise the high-risk HPVs (reviewed in 52). Amongst these E5 proteins high sequence conservation exists (reviewed in 53). Thus, the structure proposed in this study represents a large range of E5 protein structures of different genotypes of high-risk HPV-16.

ACS Paragon Plus Environment

21

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 38

Comparison with experimental data and proposed computational models

The presented model matches the experimental data of Wetherill et al. 9 in as much that the bundle fits the doughnut like shape of the protein assembly identified by TEM data of E5. TMD2 is also proposed to be pore lining and Leu-45 and -49 facing the lumen due to a docking approach applied 9

. The maximum pore radius calculated from that bundle is reported as of 5 Å. Size dependent

release of fluorescein dextrans of different stokes radii confirm a channel formation of the proteins in liposomes and a putative dimeter of the pore in the range of 1.2 – 2.8 nm. Thus, the minimum radii for the bundles in this study match rather well with the experimentally suggested data by Wetherill et al..

The present study suggests weak ion selectivity of the protein. Up to now many of the viral channel forming proteins are experimentally found to be not very selective. The overall architecture of ion channels, that hydrophilic residues of each monomer align so that rings of residues of the same kind are formed around the central pore axis, is similar to those of the ligand gated ion channels 48. These channels are known for their high selectivity for the activating ligand, but low selectivity for the ions enabling to permeate through the pore.

Mechanistic features

At a visual level, the proposed bundle model by Wetherill et al. is that of an assembly of more or less straight helices which resembles that of bundle-P4b. Due to correlation of the structure with

ACS Paragon Plus Environment

22

Page 23 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

those proposed in this study, the bundle should rather be a closed one and comparable with bundleP4b. Since helices of bundle-P1a are tilted, a tilt motion is suggested as a gating mechanism. The open-like structures used for characterizing E5 in this study reflect interim structures. They possibly exist temporarily and are screened in vivo at body temperature. Data from FCA allow to propose that high flexibility and large scale dynamics of the bundle structure support the continuous existence of a water column within the pore. It may be speculated further that the dynamics may also be a prerequisite for the bundle to ‘close’. Thus, gating could be entropically driven.

The polytopic membrane protein p7 of HCV with two TMDs shows similar behavior of asymmetric dynamics of the individual TDMs in either the monomer or the hexameric bundle 44. It is still up to further investigations in as much gating involves any concerted dynamics at all. Closing of the bundle structure is considered to be slower than the 1 ns. Thus the bundle structures used to perform PMF calculations reflect reliable open-like structures.

Selectivity

Application of PMF calculations on e.g. gramicidin A channel 54 and simulations of ion permeating through potassium channels under applied voltages

55

have shown to give reliable results. PMF

calculations have been also been performed with known structures of ion channels, e.g. a putative open structure of GLIC (PDB entry ID: 3EHZ reviewed in

22

48

) and compared with other viral channels (56,

). The PMF values calculated in this study are those from an open channel since

during the MD simulation the hexameric bundle developed from a ‘closed’ bundle into an ‘open’ one which is able to conduct water molecules. The values are within the same range (e.g. -0.7 to 2.8

ACS Paragon Plus Environment

23

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 38

kcal/mol) as found for similar calculations with the putative pore structure of another VCP 8a of SARS-CoV even though the sequence is not similar to the one used in this study 56. PMF curves for e.g. Na- or K-ions through the pentameric 8a-bundle are very much identical. Thus, as long as the overall alternating pattern of a hydrophilic (mostly due to serines and/or threonines) and hydrophobic stretches within the lumen of the pore is present, the calculations seem to be tolerable to different amino acids and turn out to be quite similar for many ions at comparable pore radii. It is rather then the pore radius which imposes considerably different values for different ions.

The permeation of the ions through the relatively large pore formed by E5 does not show any knock-on effect (as reported for potassium channel in

55

) but rather a single filed passing with

positive and negative ions either passing individually or in some occasions simultaneously through the pore. In this study Cl-ions pass more readily through the pore in the voltage simulations which parallels to some extent similar type of simulations with the viral channel forming protein p7 of HCV

57

. With increasing pore size selectivity of E5 vanishes. Thus, the term the proteins form a

‘pore’ is justified. It seems the virus is not in need of having channels which perform on a welldefined basis.

Conclusions

A fully in silico analysis of the viral channel forming protein E5 is presented in this study. The instrumentation of doing so is based on SSPPs, a docking approach and extended MD simulations in a most realistic biochemical environment. The structural data are analyzed using physical chemical tools such as PMF calculations and statistical methods such as FCA.

ACS Paragon Plus Environment

24

Page 25 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

All protocols to generate structural and functional data aim to lean on experimental data as they are available in the literature and contribute with novel information. Experimentally it has been found that E5 conducts ions in its hexameric form. At this stage, SSPPs are able to suggest residues which are supposed to be within a transmembrane stretch of the protein. The programs also allow to restrict the transmembrane motif to be helical. At this stage the structural features are experimentally not defined, but can be proposed by docking approaches. Furthermore, experiments show that the protein releases fluorescent dyes from liposomes and thus is proposed to join the class of VCPs. As a novel contribution to the question of how the protein is working, selectivity based on reliable structural models of the protein in its bundle architecture is given by computational analysis. It is proposed that (i) the protein forms channels selective for monovalent cations, (ii) that selectivity is assumed not to be high and (iii) the bundle could show weak rectification.

Comparison of these data with those, experimentally or computationally, for other suggest that VCPs in general seem to be low-precision, multi-functional tools.

Supporting Information Available The supporting information contains a table of the minimum pore radii for bundle-P1a under applied voltages. Figures are given with information about from MD simulations from bundle-P1a, -P1b and -P4b about (i) energies of the docked hexameric bundles as a function of distance, rotational and tilt angles, as well as RMSD values for the simulation of these bundles, (ii) calculated pore radii along the z-axis of the three bundles from simulations without voltage and those with either positive or negative voltages, (iii) visualization of the movements of the individual helices of

ACS Paragon Plus Environment

25

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 38

the three bundles based on FCA, (iv) projected displacement of the first two FCA modes of the three bundles over the entire duration of the MD simulation as well as over the last 20 ns.

Acknowledgment

WBF thanks the Ministry of Science and Technology, Taiwan, (NSC-101-2112-M-010-002-MY3) for financial support. DRM acknowledges a Taiwan scholarship from the Ministry of Education of Taiwan. We are grateful to the National Center for High-performance Computing for computer time and facilities.

References

(1)

Aldrich, R. W. A new standard: a review of Handbook of Ion Channels. J. Gen. Physiol.

2015, 146, 119-121. (2)

Roux, B.; Allen, T.; Berneche, S.; Im, W. Theoretical and computational models of

biological ion channels. Quart. Rev. Biophys. 2004, 37, 15-103. (3)

Maffeo, C.; Bhattacharya, S.; Yoo, J.; Wells, D.; Aksimentiev, A. Modeling and simulation

of ion channels. Chem. Rev. 2012, 112, 6250-6284. (4)

Berti, C.; Furini, S.; Gillespie, D.; Boda, D.; Eisenberg, R. S.; Sangiorgi, E.; Fiegna, C.

Three-dimensional Brownian dynamics simulator for the study of ion permeation through membrane pores. J. Chem. Theory Comput. 2014, 10, 2911-2926.

ACS Paragon Plus Environment

26

Page 27 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(5)

The Journal of Physical Chemistry

Trick, J. L.; Aryal, P.; Tucker, S. J.; Sansom, M. S. Molecular simulation studies of

hydrophobic gating in nanopores and ion channels. Biochem. Soc. Trans. 2015, 43, 146-150. (6)

Nieva, J. L.; Madan, V.; Carrasco, L. Viroporins: structure and biological functions. Nat.

Rev. Microbiol. 2012, 10, 563-574. (7)

Fischer, W. B.; Li, L.-H.; Mahato, D. R.; Wang, Y.-T.; Chen, C. Viral channel proteins in

intracellular protein-protein communication: Vpu of HIV-1, E5 of HPV16 and p7 of HCV. Biochim. Biophys. Acta 2014, 1838, 1113-1121. (8)

Bubb, V. J.; McCance, D. J.; Schlegel, R. DNA sequence of the HPV-16 E5 ORF and the

structural conservation of its encoded protein. Virology 1988, 163, 243-246. (9)

Wetherill, L. F.; Holmes, K. K.; Verow, M.; Müller, M.; Howell, G.; Harris, M.; Fishwick,

C.; Stonehouse, N.; Foster, R.; Blair, G. E., et al. High-risk human papillomavirus E5 oncoprotein displays channel-forming activity sensitive to small-molecule inhibitors. J. Virol. 2012, 86, 53415351. (10)

zur Hausen, H. Papillomaviruses in anogenital cancer as a model to understand the role of

viruses in human cancers. Cancer Res. 1989, 49, 4677-4681. (11)

zur Hausen, H. Papillomaviruses and cancer: from basic studies to clinical application. Nat.

Rev. Cancer 2002, 2, 242-350. (12)

Conrad, M.; Bubb, V. J.; Schlegel, R. The human papillomavirus type 6 and 16 E5 proteins

are membrane-associated proteins which associate with the 16-kilodalton pore-forming protein. J. Virol. 1993, 67, 6170-6178. (13)

Krawczyk, E.; Suprynowicz, F. A.; Sudarshan, S. R.; Schlegel, H. B. Membrane orientation

of the human papillomavirus type 16 E5 oncoprotein. J. Virol. 2010, 84, 1696-1703.

ACS Paragon Plus Environment

27

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(14)

Page 28 of 38

Gieswein, C. E.; Sharom, F. J.; Wildeman, A. G. Oligomerization of the E5 protein of

human papillomavirus type 16 occurs through multiple hydrophobic regions. Virology 2003, 313, 415-426. (15)

Venuti, A.; Paolini, F.; Nasir, L.; Corteggio, A.; Roperto, S.; Campo, M. S.; Borzacchiello,

G. Papillomavirus E5: the smallest oncoprotein with many functions. Mol. Cancer 2011, 10, 140. (16)

Ullman, C. G.; Haris, P. I.; Kell, B.; Cason, J.; Jewers, R. J.; Best, J. M.; Emery, V. C.;

Perkins, S. J. Hypothetical structure of the membrane-associated E5 oncoprotein of human papillomavirus type 16. Biochem. Soc. Trans. 1994, 22, 439S. (17)

Nath, R.; Mant, C. A.; Kell, B.; Cason, J.; Bible, J. M. Analysis of variant human

papillomavirus type-16 E5 proteins for their ability to induce mitogenesis of murine fibroblasts. Cancer Cell Int. 2006, 6, 19. (18)

Kell, B.; Jewers, R. J.; Cason, J.; Pakarian, F.; Kaye, J. N.; Best, J. M. Detection of E5

oncoprotein in human papillomavirus type 16-positive cervical scrapes using antibodies raised to synthetic peptides. J. Gen. Virol. 1994, 75, 2451-2456. (19)

Alonso, A.; Reed, J. Modelling of the human papillomavirus type 16 E5 protein. Biochim.

Biophys. Acta 2002, 1601, 9-18. (20)

Bauer-Hofmann, R.; Borghouts, C.; Auvinen, E.; Bourda, E.; Rösl, F.; Alonso, A. Genomic

cloning and characterization of the nonoccupied allele corresponding to the integration site of human papillomavirus type 16 DNA in the cervical cancer cell line SiHa. Virology 1996, 217, 3341. (21)

Stoler, M. H.; Rhodes, C. R.; Whitbeck, A.; Wolinsky, S. M.; Chow, L. T.; Broker, T. R.

Human papillomavirus type 16 and 18 gene expression in cervical neoplasias. Human Pathology 1992, 23, 117-127.

ACS Paragon Plus Environment

28

Page 29 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(22)

The Journal of Physical Chemistry

Fischer, W. B.; Kalita, M. M.; Heermann, D. Viral Channel forming proteins - how to

assemble and depolarize lipid membranes in silico. Biochim. Biophys. Acta 2016, 1858, 1710-1721. (23)

Lange, O. F.; Grubmüller, H. Full correlation analysis of conformational protein dynamics.

Proteins 2008, 70, 1294-1312. (24)

Flores, E. R.; Allen-Hoffmann, B. L.; Lee, D.; Sattler, C. A.; Lambert, P. F. Establishment

of the human papillomavirus type 16 (HPV-16) life cycle in an immortalized human foreskin keratinocyte cell line. Virology 1999, 262, 344-354. (25)

Hwang, E.-S.; Nottoli, T.; Dimaio, D. The HPV16 E5 protein: expression, detection, and

stable complex formation with transmembrane proteins in COS cells. Virology 1995, 211, 227-233. (26)

Yang, D. H.; Wildeman, A. G.; Sharom, F. J. Overexpression, purification, and structural

analysis of the hydrophobic E5 protein from human papillomavirus type 16. Protein Expr. Purif. 2003, 30, 1-10. (27)

Krogh, A.; Larsson, B.; von Heijne, G.; Sonnhammer, E. L. L. Predicting transmembrane

protein topology with a hidden Markov model: application to complete genomes. J. Mol. Biol. 2001, 305, 567-580. (28)

Cserzö, M.; Wallin, E.; Simon, I.; Von Heijne, G.; Elofsson, A. Prediction of

transmembrane alpha-helices in procariotic membrane proteins: the Dense Alignment Surface method. Protein Eng. 1997, 10, 673-676. (29)

Jones, D. T. Protein secondary structure prediction based on position-specific scoring

matrices. J. Mol. Biol. 1999, 292, 195-202. (30)

Cao, B.; Porollo, A.; Adamczak, R.; Jarrell, M.; Meller, J. Enhanced recognition of protein

transmembrane domains with prediction-based structural profiles. Bioinformatics 2006, 22, 303309.

ACS Paragon Plus Environment

29

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(31)

Page 30 of 38

Käll, L.; Krogh, A.; Sonnhammer, E. L. A combined transmembrane topology and signal

peptide prediction method. J. Mol. Biol. 2004, 338, 1027-1036. (32)

Viklund, H.; Elofsson, A. OCTOPUS: improving topology prediction by two-track ANN-

based preference scores and an extended topological grammar. Bioinformatics 2008, 24, 16621668. (33)

Kahsay, R. Y.; Gao, G.; Liao, L. An improved hidden Markov model for transmembrane

protein detection and topology prediction and its applications to complete genomes. Bioinformatics 2005, 21, 1853-1858. (34)

Tusnády, G. E.; Simon, I. The HMMTOP transmembrane topology prediction server.

Bioinformatics 2001, 17, 849-850. (35)

Cheng, J.; Baldi, P. Three-stage prediction of protein β-sheets by neural networks,

alignments and graph algorithms. Bioinformatics 2005, 21, Suppl. 1, i75-i84. (36)

Subramania, A.; Floudas, C. A. β-sheet topology prediction with high precision and recall

for β and mixed α/β proteins. PLoS One 2012, 7, e32461. (37)

Soto, C. S.; Fasnacht, M.; Zhu, J.; Forrest, L.; Honig, B. Loop modeling: sampling, filtering,

and scoring. Proteins 2008, 70, 834-843. (38)

Krüger, J.; Fischer, W. B. Exploring the conformational space of Vpu from HIV-1: a

versatile adaptable protein. J. Comput. Chem. 2008, 29, 2416-2424. (39)

Krüger, J.; Fischer, W. B. Assembly of viral membrane proteins. J. Chem. Theory Comput.

2009, 5, 2503-2513. (40)

Hsu, H.-J.; Fischer, W. B. In silico investigations of possible routes of assembly of ORF 3a

from SARS-CoV. J. Mol. Mod. 2011, 18, 501-514.

ACS Paragon Plus Environment

30

Page 31 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(41)

The Journal of Physical Chemistry

Adamczak, R.; Porollo, A.; Meller, J. Accurate prediction of solvent accessibility using

neural networks based regression. Proteins 2004, 56. (42)

Drozdetskiy, A.; Cole, C.; Procter, J.; Barton, G. J. JPred4: a protein secondary structure

prediction server. Nucleic. Acids Res. 2015, 43, W389-W394. (43)

Chandrasekhar, I.; Kastenholz, M.; Lins, R. D.; Oostenbrink, C.; Schuler, L. D.; Tieleman,

D. P.; van Gunsteren, W. F. A consistent potential energy parameter set for lipids: dipalmitoylphosphatidylcholine as a benchmark of the GROMOS96 45A3 force field. Eur. Biophys. J. 2003, 32, 67-77. (44)

Kalita, M. M.; Fischer, W. B. Asymmetric dynamics of ion channel forming proteins -

hepatitis C virus (HCV) p7 bundles. Biochim Biophys Acta 2016, 1858, 1462-1470. (45)

Smart, O. S.; Neduvelil, J. G.; Wang, X.; Wallace, B. A.; Sansom, M. S. P. Hole: A program

for the analysis of the pore dimensions of ion channel structural models. J. Mol. Graph. 1996, 14, 354 - 360. (46)

Nishimura, K.; Kim, S.; Zhang, L.; Cross, T. A. The closed state of a H+ channel helical

bundle combining precise orientational and distance restraints from solid state NMR. Biochemistry 2002, 41, 13170-13177. (47)

Hu, J.; Asbury, T.; Achuthan, S.; Bertram, R.; Quine, J. R.; Fu, R.; Cross, T. A. Backbone

structure of the amantadine-blocked trans-membrane domain M2 protein channel from influenza A virus. Biophys. J. 2007, 92, 4335-4343. (48)

Hilf, R. J. C.; Dutzler, R. Structure of a potentially open state of a proton-activated

pentameric ligand-gated ion channel. Nature 2009, 457, 115-119.

ACS Paragon Plus Environment

31

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(49)

Page 32 of 38

Li, L.-H.; Hsu, H.-J.; Fischer, W. B. Qualitative computational bioanalytics: Assembly of

viral channel-forming peptides around mono and divalent ions. Biochem. Biophys. Res. Commun. 2013, 442, 85-91. (50)

Park, S. H.; Mrse, A. A.; Nevzorov, A. A.; Mesleh, M. F.; Oblatt-Montal, M.; Montal, M.;

Opella, S. J. Three-dimensional structure of the channel-forming trans-membrane domain of virus protein "u" (Vpu) from HIV-1. J. Mol. Biol. 2003, 333, 409-424. (51)

Bravo, I. G.; Alonso, A. Mucosal human papillomaviruses encode four different E5 proteins

whose chemistry and phylogeny correlate with malignant or benign growth. J. Virol. 2004, 78, 13613-13626. (52)

de Villiers, E. M.; Fauquet, C.; Broker, T. R.; Bernard, H. U.; zur Hausen, H. Classification

of papillomaviruses. Virology 2004, 324, 17-27. (53)

DiMaio, D.; Petti, L. M. The E5 proteins. Virology 2013, 445, 99-114.

(54)

Allen, T. W.; Andersen, O. S.; Roux, B. Molecular dynamics - potential of mean force

calculations as a tool for understanding ion permeation and selectivity in narrow channels. Biophys. Chem. 2006, 124, 251-267. (55)

Jensen, M. Ø.; Jogini, V.; Eastwood, M. P.; Shaw, D. E. Atomic-level simulation of current-

voltage relationships in single-file ion channels. J. Gen. Physiol. 2013, 141, 619-632. (56)

Hsu, H.-J.; Lin, M.-H.; Schindler, C.; Fischer, W. B. Structure based computational

assessment of channel properties of assembled ORF-8a from SARS-CoV. Proteins 2015, 83, 300308. (57)

Wang, Y.-T.; Schilling, R.; Fink, R. H. A.; Fischer, W. B. Ion-dynamics in hepatitis C virus

p7 helical transmembrane domains - a molecular dynamics simulation study. Biophys. Chem. 2014, 192, 33-40.

ACS Paragon Plus Environment

32

Page 33 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Fig. 1: (A) Identification of helical motifs in the sequence of E5 protein using a series of secondary structure prediction programs for α-helical motifs with a TMD, TMHMM, DASTM, PSIPRED, MINNOU, PHOBIUS, TOPO2, OCTOPUS, TMMOD and HMMTOP. Predicted residues to adopt a helical motif are marked with apteryx and an orange bar. Programs BETApro and BeST are used for predicting β-fold. (B) Structural model according to P1a (see Table 1) with TMD2 facing the pore. Ideal helices TMD1 and TMD2 are docked first. TMD3 is then docked onto the first two TMDs to generate the monomer. The monomer is then copied six times and the bundle formed by a concerted docking approach. Hydrophilic residues Ser-35, -37, -41 and -52 (red) as well as Thr-38 and -49 (yellow) are highlighted. All other residues of TMD2 are highlighted in orange. Atoms are depicted in van der Waals representation and two monomers are removed to view inside the lumen of the pore. 193x216mm (300 x 300 DPI)

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Fig. 2: Hexameric bundles derived from three protocols. Bundle-P1a, bundle-P1b and bundle-P4b, are shown in a side view at the beginning (0 ns) and the end of the MD simulation (100 ns). Protein backbones are shown in yellow ribbon representation and hydrophilic residues such as serines and threonines highlighted in red stick modus. The van der Waals surface of the entire protein is shown in light grey. Water molecules are shown in blue van der Waals surface representation. 181x118mm (300 x 300 DPI)

ACS Paragon Plus Environment

Page 34 of 38

Page 35 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Fig. 3: Visualization of the first eigenvector of the structures calculated from a full correlation analysis (FCA), FCA1, of bundle-P1a (top row) over the entire duration of the simulations (100 ns) and the last 20 ns (80 – 100 ns) (first and last image, respectively). The color of the bundles follows the coding for B-factor calculations with blue represents low B-factor and red represents a high value. Averaged B-factors over all Cα-atoms and (i) averaged over 100 ns are: bundles-P1a ((37.3 ±11.7) Å2), bundle-P1b ((31.8 ± 5.8) Å2), bundle-P4b ((25.2 ± 7.9) Å2), as well as (ii) averaged over the last 20 ns: bundle-P1a (18.7 ±4.2 Å2), bundle-P1b (17.1 ± 5.8 Å2), bundle-P4a (14.7 ± 2.0 Å2). Projection of the trajectory on the first FCA mode with the second FCA mode of bundle-P1a (MI = 0.35) from the entire duration of the simulation as well as the last 20 ns (MI = 0.05) (first and second plot, respectively of top row). Similar images and plots for bundle-P1b (MI = 0.44 (100 ns), MI = 0.01 (last 20 ns) (middle row) as well as bundle-P4b (MI = 0.54 (100 ns), MI = 0.00 (last 20 ns) bottom row). 190x158mm (300 x 300 DPI)

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Fig. 4: Potential of Mean Force (PMF) calculation for Ca- (blue), Na- (red), K- (green), and Cl-ions (black) along the pore axis of bundles-P1a and -P4b. The frames for PMF calculation were taken from 20ns simulation while both the hexamers having water molecules within the entire length of the pore. The backbone of the bundles is shown in ribbon mode (white) and hydrophilic residues of TMD2 (yellow) lining the pore of hexamer are highlighted in stick mode (red). The PMF values are evaluated using bootstrap analysis (No. of bins = 200). 191x154mm (300 x 300 DPI)

ACS Paragon Plus Environment

Page 36 of 38

Page 37 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Fig. 5: Number of ions permeating tough bundle-P1a during the MD simulation under applied voltages such as +/-5 mV, +/- 11 mV, +/-33 mV, +/55 mV, +/-77 mV and 99 mV. Bundle-P1a is embedded in a 1 M solution of NaCl (109 Na-ions / 127 Cl-ions). Positive voltage leads to the movement of the Na-ion from N→C, negative voltage leads to the movement of the Cl-ions from N→C. The view into the pore of the bundles at voltages +/- 5 mV, +/- 55 mV and +/- 99 mV are shown underneath the curve. The structures are the last frame of the 200 ns MD simulation. 194x232mm (300 x 300 DPI)

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Table of Content Graphics 180x131mm (300 x 300 DPI)

ACS Paragon Plus Environment

Page 38 of 38