Supplementary material
Structure / Starting / Final unrestrained (35 ns) / Unrestrained after the 3rd annealing (35 ns in total)Total / Total / Out / In / Total / Out / In
KcsA_KcsA / 25 / 23 / 2 / 0 / 21 / 5 / 1
KcsA_KirBac / 38 / 19 / 28 / 9 / 11 / 34 / 7
KcsA_MlotiK / 4 / 5 / 0 / 1 / 5 / 1 / 2
KcsA_NaK / 34 / 9 / 34 / 9 / 17 / 29 / 12
KirBac_KcsA / 25 / 2 / 23 / 0 / 10 / 18 / 3
KirBac_KirBac / 38 / 14 / 36 / 12 / 13 / 38 / 13
KirBac_MlotiK / 4 / 0 / 4 / 0 / 0 / 4 / 0
KirBac_NaK / 34 / 33 / 34 / 33 / 29 / 34 / 29
MlotiK_KcsA / 25 / 13 / 14 / 2 / 14 / 12 / 1
MlotiK_KirBac / 38 / 27 / 37 / 26 / 18 / 37 / 17
MlotiK_MlotiK / 4 / 0 / 4 / 0 / 0 / 4 / 0
MlotiK_NaK / 34 / 29 / 30 / 25 / 33 / 32 / 31
NaK_KcsA / 25 / 17 / 12 / 4 / 20 / 5 / 0
NaK_KirBac / 38 / 17 / 27 / 6 / 34 / 10 / 6
NaK_MlotiK / 4 / 15 / 4 / 15 / 15 / 4 / 15
NaK_NaK / 34 / 25 / 28 / 19 / 21 / 27 / 14
Table SI. Exchange of the water molecules between the central cavity of the channel and the bulk water. “Out” and “In” columns indicate number of the water molecules that have left and entered the cavity respectively during the 35 ns simulation in total.
Structure / Residues / # atoms / # watermolecules / # lipid
molecules / # anions/cations / cell size, Å
KcsA / Ser22-Gly116 / 69470 / 11000 / 229 / 46(Cl-)/38(K+) / 97x97x78
KirBac / Thr43-Arg137 / 68551 / 11065 / 220 / 42(Cl-)/42(K+) / 97x97x78
NaK / Lys22-Gln103 / 68279 / 11066 / 222 / 39(Cl-)/41(Na+)+1(Ca2+) / 97x97x78
MlotiK / Ala126-Val195 / 68049 / 10979 / 220 / 40(Cl-)/44(K+) / 97x97x78
Table SII: Description of the systems used during simulations.
Fig. S1: (A) Deviation from symmetry and the value of the symmetry restraint spring constant vs. time; (B) deviation from symmetry vs. spring constant.
Fig. S2: RMSD of the entire trace calculated against the initial homology models, Ca atoms used for alignment and for RMSD calculation
Fig. S3: RMSD against target x-ray structure, all heavy atoms used for alignment and for RMSD calculation.
Fig. S4. RMSD calculation for individual monomers, using the CA atoms for alignment and calculation. (A) chain A; (B) chain B; (C) chain C; (D) chain D.
(A)
(B)
(C)
(D)
Fig. S5: Presence of ions in the selectivity filter
Fig. S6: Correlation between HMs accuracy, measured as RMSD against x-ray structure, and stability, measured as RMSD against initial frame of each trajectory, during different stages of unrestrained MD. (A) after the third annealing; (B) after the second annealing; (C) after the first annealing; (D) during the last 8ns of unrestrained simulation. Numbers correspond to different HMs:
2 / NaK_MlotiK / 5 / NaK_KcsA / 8 / KcsA_NaK / 11 / NaK_KirBac
3 / KirBac_MlotiK / 6 / KirBac_KcsA / 9 / KirBac_NaK / 12 / KcsA_KirBac
Fig. S7: RMSD of all Ca atoms as a function of simulation time (ns), calculated using the x-ray structure of the target channel as reference, for x-ray structures and HMs of each protein: MlotiK (first row), KcsA second row), NaK (third row), and KirBac (fourth row). The color scheme is the same as in Fig. 3
Sample NAMD configuration file for the symmetry annealing simulation steps:
set topologyFilePath {/data/amilac/0symm_andriy/KcsA_KcsA} ; # Folder with topology/parameters files
set structureFileNameBase {KcsA_KcsA_ornt_popc_wat_fi_salt}
set startupFileNameBase {sim-600_KcsA_KcsA_ornt_popc_wat_fi_salt__no-cmap_6645000_alnCmass} ; # Starting coordinates filename
set outputFileNameBase {symm-300_KcsA_KcsA_ornt_popc_wat_fi_salt__cmap} ; # Output filenames will be constructed with this filename base plus extension
set firstTimestep 0 ; # The very first timestep
set lastTimestep 1000000 ; # The very last timestep
set workingFolder {/data/amilac/0symm_andriy/KcsA_KcsA} ; # The folder for all the simulation files
set indicesFilePathName {/data/amilac/0symm_andriy/KcsA_KcsA/KcsA_KcsA_ornt_popc_wat_fi_salt_indices.txt}; # File with indices of the atoms that should be symmetrized during the annealing
set targetUpdateSteps 1000
set monomers 4; # Number of monomers
set rotation_axis z; #Rotation symmetry axis - x,y or z
set rotation_direction 1; # 1: Clockwise; -1: CounterClockwise (looking in the positive z direction)
set forceUpdateSteps 1; # The number of timesteps between the force direction/amount update. If more than 1, the force is applied only every i-th step and scaled up $forceUpdateSteps times while no force is applied in between.
set minimalConstant 0.001; # Initial/minimal value of harmonic spring constant
set maximalConstant 10.0; # maximal value of harmonic spring constant
set multiplierConstant 1.1; # The multiplier for exponential growth/decay of the spring constant
set targetRmsd 0.01; # The annealing will increase/decrease the force constant starting with the $initialConstant and $multiplierConstant multiplier to provide linear decrease of RMSD from the starting to the target value
set linearRmsdDecrease 1; # If 1 - the spring coefficient will grow/decrease exponentially trying to decrease/grow RMSD from the initial to the target value in time. If 0 - the spring coefficient will grow exponentially until the target RMSD is reached, then coefficient will grow/decrease to maintain the target RMSD.
# ______Procedures declaration
proc saveXYZ {xyzList fileName {title ""}} {
# Procedure receives a list of the coordinate triplets and saves them as a file with fileName in crd-like format
set fID [open $fileName w]
puts $fID "TITLE : $title"
foreach {x y z} [join [join $xyzList]] {puts $fID "$x $y $z"}
close $fID
}
proc cartesian2polar {xy} {
# Procedure converts array (length 2) from xy cartesian coordinates to rt polar coordinates
# It can also convert a list of such pairs {{xy} {xy} ... {xy}}
# If instead of pair(s) of coordinates it gets them by 3, it considers tham as {x y z}, coverts x and y to polar and leaves z as it was
if {[llength [lindex $xy 0]] > 1} {
# This is a list of pairs
set rtList {}
set xyList $xy
# Go through the list and convert coordinates
if {[llength [lindex $xyList 0]]==2} {
# 2D coordinates
foreach xy $xyList {
set rt {}
lappend rt [expr {pow((pow([lindex $xy 0],2)+pow([lindex $xy 1],2)),0.5)}]
lappend rt [expr {atan2([lindex $xy 1],[lindex $xy 0])}]
lappend rtList $rt
}
} else {
# 3D coordinates
foreach xy $xyList {
set rt {}
lappend rt [expr {pow((pow([lindex $xy 0],2)+pow([lindex $xy 1],2)),0.5)}]
lappend rt [expr {atan2([lindex $xy 1],[lindex $xy 0])}]
lappend rt [lindex $xy 2]
lappend rtList $rt
}
}
return $rtList
} else {
# This is a single pair
# Go through the list and convert coordinates
if {[llength $xy]==2} {
# 2D coordinates
set rt {}
lappend rt [expr {pow((pow([lindex $xy 0],2)+pow([lindex $xy 1],2)),0.5)}]
lappend rt [expr {atan2([lindex $xy 1],[lindex $xy 0])}]
} else {
# 3D coordinates
set rt {}
lappend rt [expr {pow((pow([lindex $xy 0],2)+pow([lindex $xy 1],2)),0.5)}]
lappend rt [expr {atan2([lindex $xy 1],[lindex $xy 0])}]
lappend rt [lindex $xy 2]
}
return $rt
}
}
proc polar2cartesian {rt} {
# Procedure converts array (length 2) from rt polar coordinates to xy cartesian coordinates
# It can also convert a list of such pairs {{rt} {rt} ... {rt}}
# If instead of pair(s) of coordinates it gets them by 3, it considers tham as {r t z}, coverts r and t to cartesian and leaves z as it was
if {[llength [lindex $rt 0]] > 1} {
# This is a list of pairs
set xyList {}
set rtList $rt
# Go through the list and convert coordinates
if {[llength [lindex $rtList 0]]==2} {
# 2D coordinates
foreach rt $rtList {
set xy {}
lappend xy [expr {[lindex $rt 0]*cos([lindex $rt 1])}]
lappend xy [expr {[lindex $rt 0]*sin([lindex $rt 1])}]
lappend xyList $xy
}
} else {
# 3D coordinates
foreach rt $rtList {
set xy {}
lappend xy [expr {[lindex $rt 0]*cos([lindex $rt 1])}]
lappend xy [expr {[lindex $rt 0]*sin([lindex $rt 1])}]
lappend xy [lindex $rt 2]
lappend xyList $xy
}
}
return $xyList
} else {
# This is a single pair
if {[llength [lindex $rt 0]]==2} {
# 2D coordinates
set xy {}
lappend xy [expr {[lindex $rt 0]*cos([lindex $rt 1])}]
lappend xy [expr {[lindex $rt 0]*sin([lindex $rt 1])}]
lappend xyList $xy
} else {
# 3D coordinates
set xy {}
lappend xy [expr {[lindex $rt 0]*cos([lindex $rt 1])}]
lappend xy [expr {[lindex $rt 0]*sin([lindex $rt 1])}]
lappend xy [lindex $rt 2]
}
return $xy
}
}
proc symmRmsdCoords {coordsStart monomers rotation_axis rotation_direction} {
set pi [expr {acos(-1)}]
set monomerAtoms [expr {int([llength $coordsStart]/$monomers)}]
# Rotate all the monomers
set coordsStartPolar [cartesian2polar $coordsStart]
set coordsRot [lrange $coordsStartPolar 0 [expr {$monomerAtoms-1}]]
for {set m 1} {$m < $monomers} {incr m} {; # Go through the monomers and rotate them
set rotationAngle [expr {-2.0*$pi*$m*$rotation_direction/$monomers}]; # Each monomer will be rotated to overlap with the 0'th one
set monomerAtomsStart [expr {int($m*$monomerAtoms)}]
set monomerAtomsEnd [expr {int(($m+1)*$monomerAtoms)}]
for {set ma $monomerAtomsStart} {$ma<$monomerAtomsEnd} {incr ma} {; # Go through the atoms of the current monomer and rotate them
foreach {x y z} [lindex $coordsStartPolar $ma] {}
lappend coordsRot "$x [expr {$y+$rotationAngle}] $z"
}
}
set coordsRot [polar2cartesian $coordsRot]
# Average all the monomers
set coordsAvg [lrange $coordsRot 0 [expr {$monomerAtoms-1}]]
for {set m 1} {$m < $monomers} {incr m} {; # Go through the monomers and add coordinates
set monomerAtomsStart [expr {int($m*$monomerAtoms)}]
set monomerAtomsEnd [expr {int(($m+1)*$monomerAtoms)}]
set coordsAvgNew {}
for {set ma 0} {$ma<$monomerAtoms} {incr ma} {; # Go through the atoms of the current monomer and add them
foreach {x0 y0 z0} [lindex $coordsAvg $ma] {}
foreach {x y z} [lindex $coordsRot [expr {$ma+int($m*$monomerAtoms)}]] {}
lappend coordsAvgNew "[expr {$x0+$x}] [expr {$y0+$y}] [expr {$z0+$z}]"
}
set coordsAvg $coordsAvgNew
}
set coordsAvgNew {}
for {set ma 0} {$ma<$monomerAtoms} {incr ma} {; # Go through the atoms of the current monomer and divide them by the number of monomers
foreach {x y z} [lindex $coordsAvg $ma] {}
lappend coordsAvgNew "[expr {double($x/$monomers)}] [expr {double($y/$monomers)}] [expr {double($z/$monomers)}]"
}
set coordsAvg [list $coordsAvgNew]
for {set m 1} {$m < $monomers} {incr m} {; # Go through the monomers and duplicate monomer coordinates
lappend coordsAvg $coordsAvgNew
}
set coordsAvg [join $coordsAvg]
# Compare deviation of all the monomers from the average
set symmRmsd 0
foreach xyzRor $coordsRot xyzAvg $coordsAvg {
foreach {xRot yRot zRot} $xyzRor {xAvg yAvg zAvg} $xyzAvg {}
set symmRmsd [expr {$symmRmsd+($xRot-$xAvg)*($xRot-$xAvg)+($yRot-$yAvg)*($yRot-$yAvg)+($zRot-$zAvg)*($zRot-$zAvg)}]
}
set symmRmsd [expr {pow($symmRmsd/[llength $coordsRot],0.5)}]
# Duplicate and rotate back the first monomer
set coordsAvg [lrange $coordsAvg 0 [expr {$monomerAtoms-1}]]
set coordsAvg [cartesian2polar $coordsAvg]
set coordsAvgNew $coordsAvg
for {set m 1} {$m < $monomers} {incr m} {; # Go through the monomers and rotate them
set rotationAngle [expr {2.0*$pi*$m*$rotation_direction/$monomers}]; # Each monomer will be rotated to overlap with the 0'th one
set monomerAtomsStart [expr {int($m*$monomerAtoms)}]
set monomerAtomsEnd [expr {int(($m+1)*$monomerAtoms)}]
for {set ma 0} {$ma<$monomerAtoms} {incr ma} {; # Go through the atoms of the current monomer and rotate them
foreach {x y z} [lindex $coordsAvg $ma] {}
lappend coordsAvgNew "$x [expr {$y+$rotationAngle}] $z"
}
}
set coordsAvg [polar2cartesian $coordsAvgNew]
return [lappend symmRmsd $coordsAvg]
}
# ______Computational part
set symmRmsdStarting {}
set symmRmsdRequired $targetRmsd
set cycleConstant $minimalConstant
set symmTargetCoordinates {}
set restartFileNameBase $startupFileNameBase; # Define filename for the initial coordinates
# ______Set simulation parameters
structure $structureFileNameBase.psf
coordinates $structureFileNameBase.pdb
bincoordinates $restartFileNameBase.coor
binvelocities $restartFileNameBase.vel
extendedSystem $restartFileNameBase.xsc
parameters $topologyFilePath/par_all27_prot_lipid__CMAP.inp
paraTypeCharmm on
outputName $outputFileNameBase
XSTfile $outputFileNameBase.xst
outputEnergies 1000
outputTiming 1000
outputPressure 1000
dcdFreq 1000
restartFreq 1000
xstFreq 1000
wrapAll on
wrapNearest on
numsteps [expr {$lastTimestep-$firstTimestep}]
firstTimestep $firstTimestep
timestep 1
nonBondedFreq 2
fullElectFrequency 4
stepsPerCycle 20
switching on
switchDist 10
cutoff 12
pairlistdist 13.5
cellBasisVector1 97 0 0
cellBasisVector2 0 97 0
cellBasisVector3 0 0 78
cellOrigin 00.00 00.00 00.00
Pme on
PmeGridsizeX 100
PmeGridsizeY 100
PmeGridsizeZ 80
exclude scaled1-4
1-4scaling 1.0
langevin on
langevinDamping 1
langevinTemp 303.15
langevinHydrogen no
langevinPiston on
langevinPistonTarget 1.01325
langevinPistonPeriod 200
langevinPistonDecay 100
langevinPistonTemp 303.15
SurfaceTensionTarget 40
useGroupPressure yes
useFlexibleCell yes
DCDUnitCell yes
useConstantRatio yes
# ______Force computation
set i0 $firstTimestep
# Get the indices of the symmetrized atoms
cd $workingFolder
set il {}; set fid [open $indicesFilePathName r]; while {[gets $fid line] > 0} {lappend il $line}; close $fid; set il [join $il]
set numatoms [llength $il]
set namdIndexList {}; foreach ila $il {lappend namdIndexList [expr {$ila+1}]}; set namdIndexList [lsort -integer $namdIndexList]
tclForces on
tclForcesScript {
print "Load Forces script started"
foreach name $namdIndexList {addatom $name}
proc calcforces {} {
global namdIndexList firstTimestep targetUpdateSteps monomers rotation_axis rotation_direction minimalConstant maximalConstant multiplierConstant targetRmsd linearRmsdDecrease outputFileNameBase cycleConstant symmRmsdStarting symmRmsdCycle symmRmsdRequired symmTargetCoordinates lastTimestep forceUpdateSteps
set coordsArrayFlag 0
set step [getstep]
if {($step%$targetUpdateSteps==0)||($step==$firstTimestep)} {; # Update the symmetric target
loadcoords coordsArray
set coordsArrayFlag 1
# Estimate the RMSD from symmetric average and the symmetric target coordinates
set coordsCycle {}; foreach name $namdIndexList {lappend coordsCycle $coordsArray($name)}
foreach {symmRmsdCycle symmTargetCoordinates} [symmRmsdCoords $coordsCycle $monomers $rotation_axis $rotation_direction] {break}