MD simulations

All the inputs from molecular dynamics (MD) simulations were generated input generator tool - Membrane Builder in CHARMM-GUI (https://www.charmm-gui.org/?doc=input/membrane.bilayer). You should create an account in CHARMM-GUI to use all tools.

However, there are some preparative steps before input generation.

Preparative steps

1) Prediction of protein positioning in lipid bilayers: The structure of Kir2.1 without the ligands (K⁺ and Sr²⁺ ions) are submitted to a prediction of its positioning in lipid bilayers using PPM 2.0 web server (https://opm.phar.umich.edu/ppm_server2). The .pdb filed resulted will be used for the preparation of the system in CHARMM-GUI.

2) Evaluate the protonation of side chains: The pdb file of structure of Kir2.1 without the ligands needs to be converted in PQR using the a server (https://server.poissonboltzmann.org/pdb2pqr). Upload the PDB file and run as default. Download the .pqr and .log files. For the evaluation of protonation, it was used an in-house python script (see the end of file). To run this script: python protonationPDB2PQR.py file.log file.pqr 7 > protonation.txt. 7 is the pH that will be used for MD simulations. In the file protonation.txt there is a list of amino acids that should be protonated in the preparation of inputs in the CHARMM-GUI (Chain A/PROA: Glu138, Glu224, Glu299, His53, His110, His197, His221, His226, His271, His324, His335, His345; Chain B/PROB: Glu138, Glu293, His53, His110, His197, His221, His226, His271, His324, His335, His34; Chain C/PROC: Glu138, Lys185, Glu293, Glu299, His53, His110, His197, His221, His226, His271, His324, His335, His345; Chain D/PROD: Glu138, His53, His110, His197, His221, His226, His271, His324, His335, His345.

Preparation of inputs for MD simulations

1) Upload the PDB file that was generated using PPM server (preparative step 1) in Input Generator Tool - Membrane Builder in CHARMM-GUI (https://www.charmm-gui.org/?doc=input/membrane.bilayer) and go to next step (select model/chain).

2) Choose the 4 chains and go to the next step (Manipulate PDB);

3) Click in the box of protonation state and includes all the protonation previously predicted in the preparative step 2. After, click in the box of disulfide bonds and include PROA 122 - PROA 154 / PROB 122 - PROB 154 / PROC 122 - PROC 154 / PROD 122 - PROD 154. Go to the next step (Generate PDB and Orient Molecule);

4) In the orientation options, click in Use PDB orientation and go to the next step (Calculate cross-sectional area);

5) Click on view structure at step2_oriented.pdb to evaluate the good orientation of the membrane in relation to the protein structure. At system syze determination options insert 110 at Lenght of X and Y (initial guess). After, search for PO lipids, click in the arrow and insert 140 in the upperleflter ratio and 148 in the lowerleaft ratio for POPC lipids. Click in the show the system info. Check: Heterogenous Lipid marked, Lenght of Z based on water thickness 22.5 (Marked); Lenght of XY based on Ratios of lipid components (Marked). Go to the next step (Determine the System Size);

6) Put a concentration of 0.15 KCl (neutralizing option should be marked). Click in Calculate solvent composition. Check: Replacement Method and check lipid ring (Marked); Include Ions (Marked); Basic Ion types (Marked). Go to the next step(Build Components). If the CHARMM ends abnormally, refresh the page using F5 until it works). The following two steps are only informational, so follow to the next step (Assemble Components).

7) In input generation options, click in the box of NAMD and set temperature for 293K. Check: Force field option CHARMM36m (all the above options should be unmarked); Equilibration options: Generate grid information for PME FFT automatically (marked) and NPT ensemble (marked). Go to the next step (generate equilibration and dynamic inputs);

8) Donwload the download.tgz file to a computer/cluster. The equilibration inputs (step 6 - composed by 6 steps) and production input (step 7, composed by one step) will not run in CHARMM-GUI web server anymore. They will run in a computer/cluster containing the software NAMD. The scripts will be inside the NAMD folder after downloading and unpacking the download.tgz file

9) The equilibration inputs 6.1, 6.2, 6.3 and 6.4 were run using NAMD without any modifications. In the input 6.5, it was changed the run parameter (last parameter of the file) to 500000 (instead of 250000 as in the original file). In the input 6.6, it was changed the firsttimestep parameter to 1135000 (instead of 885000 as in the original file) and the run parameter (last parameter of the file) to 500000 (instead of 250000 as in the original file). To run namd type in the terminal namd2 +idlepoll +p6 +devices 0,1 name_of_input.inp > name_of_input..out &. +p refers to the number of processors used for this job (in this case6, choose based on what you have available on your machine). +devices refers to the GPUs used for this job (in this case GPUs 0 and 1; choose based on what you have available on your machine).

10) The production input 7, the numsteps parameter was changed to 200000000 and the run parameter was changed to 200000000 to achieve 200 ns of MD simulations.

Analysis of MD simulations

1) The analysis of Cα root mean square deviation (rmsd) was performed in the VMD. In the VMD open the step5_input.psf, open the trajectories file (step7_production.dcd) and go Extensions > Analysis > RMSD Visualizer Tool, mark the Backbone option at Atom Selection Modifiers and then click in ALIGN. Then, in atom selection write protein and segid PROA and click in RMSD and subsquently Plot result. It wil be plotted the rmsd from chain A (PROA). For the other chains, do the same procedure, changing PROA to PROB, PROC and PROD. It will open a Multiplot window and you can export the graph in the File option.

2) The analysis of Cα root mean square fluctuation (rmsf) of the residues along MD simulations was performed was perfomed using a in-house tcl script (rmsf.tcl; see the end of file). This script must be in the same folder of all files produced by MD simulations. In the VMD open the step5_input.psf, open the trajectories file (step7_production.dcd) and go Extensions > Tk Console. Then, type source rmsf.tcl The output (rmsf.dat) will appear in the folder and it is numbered 1 to 1291 (there is no separation between the chains), so it is needed to separate the respective chains: Chain A: 1 to 325: Chain B: 326-646; Chain C: 647- 967; Chain D: 968-1291.

3) The analysis of the distances between the K64 and R67 residues to the membrane (plane of the phosphate groups from the inner membrane) was performed was perfomed using a in-house tcl script (distance_membrane.tcl; see the end of file). This script must be in the same folder of all files produced by MDsimulations. In the terminal type vmd -dispdev text -e distance_membrane.tcl. Now, edit the output .dat files in terminal to get the first and and last (column 5) colums of these files (which contain the necessary information), as example: in terminal type awk '{print $1,$5}' distance_to_membrane_monA_200ns_K64.dat > K64_A.dat. In the output K64_A.dat, it will be the distance to the membrane (axis y) along the MD simulations time (axis x; in nanoseconds).

Scripts used

Python script for evaluate the protonation of the side chains (Copy and save it as protonationPDB2PQR.py file)

#!/usr/bin/env python

# -*- coding: utf-8 -*-

# Author : Rafael Borges (kindly provided by him by personal communication)

# Objetive: Read .log and .pqr from PDB2PQR output and pH desired and write

# ASP/C-/GLU which should be protonated and CYS/TYR/LYS/ARG/N+ which

# should not be protonated. For His, it says which atom is protonated.

# USAGE: Output-ProtonationPDB2PQR.py f9481qp0n3.log f9481qp0n3.pqr 7

# f9481qp0n3.log should be the log of PDB2PQR server

# f9481qp0n3.pqr should be the pqr of PDB2PQR server

# RETURNS: acid residues which should be protonated and basic residues which

# should not be protonate. For His, it says what atom is protonated.

from __future__ import print_function

import sys,os,math

log=sys.argv[1]

pqr=sys.argv[2]

pH=float(sys.argv[3])

def CalculateEuclidianDistance (v1,v2):

      if len(v1)==len(v2):

           distance = math.sqrt (sum( [(a - b) ** 2 for a, b in zip(v1, v2)] ) )

           #distance=numpy.sqrt( (list1[0]-list2[0])**2 + (list1[1]-list2[1])**2 +

(list1[2]-list2[2])**2 )

            return distance

      else:

            print ("Different vectors lengths, correct input. Exiting function

RJB_lib.CalculateEuclidianDistance.")

            exit()

with open (log) as f: frlog=f.read()

with open (pqr) as f: frpqr=f.readlines()

ireadlog=frlog.index('Group pKa model-pKa ligand atom-type')+len(' Group pKa model-pKa ligand atom-type')+1 freadlog=frlog.index('PROPKA information')

#freadlog=ireadlog+frlog[ireadlog].index('\n\n')

logsel=frlog[ireadlog:freadlog]

#print (logsel)

lres=[]

print ('Group pKa model-pKa ligand atom-type')

for l in logsel.split('\n'):

       if l=='': break

       ls=l.split()

       res,pKa=ls[0],float(ls[3])

       #print (res,pKa,pH)

       if res in ['ASP','GLU','C-'] and pKa>pH:

             print (l)

             if res=='HIS':

                  if ls[1] not in lres: lres.append(ls[1])

        if res in ['CYS','TYR','LYS','ARG','N+'] and pKa<pH: print (l)

print('\n\nAtom of HIS which should be protonated:')

dicND={}

dicNE={}

dicH={}

for l in frpqr:

     #print (l)

     if l[17:20]=='HIS':

          resn=l[21:26]

          #print (l[:-1])

          if l[13:15]=='ND':

               if resn not in dicND: dicND[resn]=[]

               dicND[resn].append((float(l[30:38]),float(l[38:46]),float(l[46:54])))

          if l[13:15]=='NE':

               if resn not in dicNE: dicNE[resn]=[]

               dicNE[resn].append((float(l[30:38]),float(l[38:46]),float(l[46:54])))

          if l[13]=='H':

               if resn not in dicH: dicH[resn]=[]

               dicH[resn].append((float(l[30:38]),float(l[38:46]),float(l[46:54])))

nd,ne=False,False

for i in range(len(dicND[resn])):

     #print (i)

     for resn in dicND:

     #print (resn)

          for H in dicH[resn]:

               if CalculateEuclidianDistance(dicND[resn][i],H)<1.5: nd=True

               else:                                                         nd=False

               if CalculateEuclidianDistance(dicNE[resn][i],H)<1.5: ne=True

               else:                                                        ne=False

               #print

(CalculateEuclidianDistance(dicND[resn][i],H),nd,CalculateEuclidianDistance(dicNE [resn][i],H),ne)

               if nd and not ne: print (resn+' ND protonated')

               if ne and not nd: print (resn+' NE protonated')

               if nd and ne:       print (resn+' ND&NE protonated')

               #print (l[13:16],CalculateEuclidianDistance(nd,h))

               #print (l[13:16],CalculateEuclidianDistance(ne,h))

tcl script for Cα root mean square fluctuation (rmsf analysis) of the residues along the MD simulations (Copy and save it as rmsf.tcl file)

set fp [open "rmsf.dat" w]

set sel [atomselect top "protein and name CA"]

for {set i 0} {$i < [$sel num]} {incr i} {

set rmsf [measure rmsf $sel first 1 last 1999 step 1]

puts $fp "[expr {$i+1}] [lindex $rmsf $i]"

}

close $fp

tcl script to measure the distances of Lys64 and Arg67 to the membrane along MD simulations (Copy and save it as distance_membrane.tcl file)

#to execute this script type vmd -dispdev text -e distance_membrane.tcl

proc calc_dist {sel1 sel2 out} {

     #Open output file to write

     set outfile [open $out w]

     set numFrame [molinfo top get numframes]

     #Defining references

     set r1 [atomselect top $sel1]

     set r2 [atomselect top $sel2]

     puts $outfile "#Time up down dist distfix"

     for {set f 0} {$f < $numFrame} {incr f} {

       #Read each frame in trajectory

       $r1 frame $f

       $r1 update

       $r2 frame $f

       $r2 update

       #Time

       #set time [format "%.2f" [expr $f * 5 * 0.05]]

       set time [expr $f /10.]

       puts $time

       #Getting Center of Mass of all PO4 (middle of the membrane)

       set p1 [measure center $r1 weight mass]

       set commemb [lindex $p1 2]

       #Defining up and down PO4 atoms

       set up []

       set down []

       set p [$r1 get {z}]

       foreach val $p {

          if { $val < $commemb } {

             lappend down $val

          } else {

             lappend up $val

         }

     }

     #Getting average of up and down PO4 atoms

     set avup [expr ([join $up +])/[llength $up]]

     set avdown [expr ([join $down +])/[llength $down]]

     #protein residues

     set pa [measure center $r2 weight mass]

     set comprot [lindex $pa 2]

     #Fixed

     set comprotfix [expr $avdown - $comprot]

     #Write to output

     puts $outfile "$time $avup $avdown $comprot $comprotfix"

    }

 close $outfile

}

#set system nopip2

#set rep ${argv}

#Load topology/structure

mol new step6.6_equilibration.pdb type pdb waitfor all

#Remove the first frame

#animate delete beg 0 end 0

#Load trajectories

mol addfile step7_production.dcd type dcd waitfor all

#Call the function

#Selections: membrane protein residues and output file

calc_dist "resname POPC and name P" "segid PROA and resid 64" "distance_to_membrane_monA_200ns_K64.dat"

calc_dist "resname POPC and name P" "segid PROA and resid 67" "distance_to_membrane_monA_200ns_R67.dat"

calc_dist "resname POPC and name P" "segid PROB and resid 64" "distance_to_membrane_monB_200ns_K64.dat"

calc_dist "resname POPC and name P" "segid PROB and resid 67" "distance_to_membrane_monB_200ns_R67.dat"

calc_dist "resname POPC and name P" "segid PROC and resid 64" "distance_to_membrane_monC_200ns_K64.dat"

calc_dist "resname POPC and name P" "segid PROC and resid 67" "distance_to_membrane_monC_200ns_R67.dat"

calc_dist "resname POPC and name P" "segid PROD and resid 64" "distance_to_membrane_monD_200ns_K64.dat"

calc_dist "resname POPC and name P" "segid PROD and resid 67" "distance_to_membrane_monD_200ns_R67.dat"

quit