Molecular Dynamics Simulations of Protein Aggregation: Protocols for Simulation Setup and Analysis with Markov State Models and Transition Networks

Abstract

Protein disorder and aggregation play significant roles in the pathogenesis of numerous neurodegenerative diseases, such as Alzheimer’s and Parkinson’s diseases. The end products of the aggregation process in these diseases are highly structured amyloid fibrils. Though in most cases, small, soluble oligomers formed during amyloid aggregation are the toxic species. A full understanding of the physicochemical forces that drive protein aggregation is thus required if one aims for the rational design of drugs targeting the formation of amyloid oligomers. Among a multitude of biophysical and biochemical techniques that are employed for studying protein aggregation, molecular dynamics (MD) simulations at the atomic level provide the highest temporal and spatial resolution of this process, capturing key steps during the formation of amyloid oligomers. Here we provide a step-by-step guide for setting up, running, and analyzing MD simulations of aggregating peptides using GROMACS. For the analysis, we provide the scripts that were developed in our lab, which allow to determine the oligomer size and inter-peptide contacts that drive the aggregation process. Moreover, we explain and provide the tools to derive Markov state models and transition networks from MD data of peptide aggregation.

Appendices

4 Appendix 1: Input Files for the MD Simulation

In the following, the five .mdp files required to follow the MD simulation protocols in Subheading 2.1 are provided.

“ ions .mdp” for adding ions

;; ions.mdp
; Run setup
integrator	= steep
emtol	= 1000
emstep	= 0.01
nsteps	= 2000
; Neighbor search
cutoff-scheme	= Verlet
pbc	= xyz
; Electrostatics and vdW
coulombtype	= PME
pme-order	= 4
fourierspacing	= 0.1
rcoulomb	= 1.2
rvdw	= 1.2

“em.mdp” for energy minimization

;; em.mdp
; Run setup
integrator	= steep
emtol	= 500
emstep	= 0.001
nsteps	= 2000
nstxout	= 100
; Neighbor searching
cutoff-scheme	= Verlet
nstlist	= 20
ns-type	= grid
pbc	= xyz
; Electrostatics
coulombtype	= PME
pme-order	= 4
fourierspacing	= 0.1
rcoulomb	= 1.2
; VdW
rvdw	= 1.2

“nvt.mdp” for first equilibration in the NVT ensemble

;; nvt.mdp
define	= -DPOSRES
; Run setup
integrator	= md
dt	= 0.002; 2 fs
nsteps	= 50000

; Output control

nstxout	= 5000
nstvout	= 5000
nstfout	= 5000
nstlog	= 500
nstenergy	= 500
; Bonds
constraints	= all-bonds
constraint-algorithm	= LINCS
lincs-order	= 4
lincs-iter	= 1
; Neighbor searching
cutoff-scheme	= Verlet
nstlist	= 20
ns-type	= grid
pbc	= xyz
; Electrostatics
coulombtype	= PME
pme-order	= 4
fourierspacing	= 0.1
rcoulomb	= 1.2
; VdW
rvdw	= 1.2
; T coupling is on
tcoupl	= v-rescale
tc-grps	= Protein Non-Protein
tau-t	= 0.1 0.1
ref-t	= 300 300
nsttcouple	= 10
; P coupling is off
pcoupl	= no
; Velocity generation
gen-vel	= yes
gen-temp	= 300
continuation	= no

“npt.mdp” for second equilibration in the NPT ensemble

;; npt.mdp
define	= -DPOSRES
; Run setup
integrator	= md
dt	= 0.002
nsteps	= 100000
; Output control
nstxout	= 5000
nstvout	= 5000
nstfout	= 5000
nstlog	= 500
nstenergy	= 500
; Bonds
constraints	= all-bonds
constraint-algorithm	= LINCS

“md.mdp” for MD production run

;; md.mdp
; Run setup
integrator	= md
dt	= 0.002
nsteps	= 500000000; 500 million
; Output control
nstxout	= 0
nstvout	= 0
nstfout	= 0
nstlog	= 2500
nstenergy	= 2500
nstxout-compressed	= 2500
compressed-x-grps	= System
; Bonds
constraints	= all-bonds
constraint-algorithm	= LINCS
lincs-order	= 4
; Neighbor searching
cutoff-scheme	= Verlet
Nstlist	= 20
ns-type	= grid
pbc	= xyz
; Electrostatics
coulombtype	= PME
pme-order	= 4
fourierspacing	= 0.1
rcoulomb	= 1.2
; VdW
rvdw	= 1.2
; T coupling is on
tcoupl	= v-rescale
tc-grps	= Protein Non-Protein
tau-t	= 0.1 0.1
ref-t	= 300 300
nsttcouple	= 10
; Pressure coupling is off
pcoupl	= Parrinello-Rahman
pcoupltype	= isotropic
tau_p	= 2.0
ref_p	= 1.0
compressibility	= 4.5e-5
; Velocity generation
gen-vel	= no
gen-temp	= 300
continuation	= yes

5 Appendix 2: Python Script for the Calculation of the Oligomerization State and Contact Maps

import numpy as np import copy as cp import subprocess as sp import os as os import shutil as sh import MDAnalysis as mdana import sys from MDAnalysis.analysis.distances import distance_array import networkx as nx import pandas as pd import mdtraj as md import matplotlib matplotlib.use("TkAgg") from matplotlib import pyplot as plt #input parameters ref_structure=sys.argv[1] traj=sys.argv[2] Min_Distance=int(sys.argv[3]) #structure parameters topology = md.load(ref_structure).topology trajectory = md.load(traj, top=ref_structure) frames=trajectory.n_frames #Number of frames chains=topology.n_chains #Number of chains atoms=int(topology.n_atoms/chains) #Number of atoms in each monomer AminoAcids = int(topology.n_residues/chains)-2 #Number of residues per chain ('-2' to not count the N- and C- cap residues as individual residues) isum=1 atoms_list=[] atomsperAminoAcid=[] residue_list=[] for residue in topology.chain(0).residues: atoms_list.append(residue.n_atoms) residue_list.append(residue) ', '.join(map(lambda x: "'" + x + "'", str(residue_list))) #The N- and C- cap residues are part of the 1st and last residue index. If no N- and C- cap residues for the protein, comment the line below using "#" del residue_list[0]; del residue_list[-1] for ii in range(len(atoms_list)): isum = isum + atoms_list[ii] atomsperAminoAcid.append(isum) atomsperAminoAcid.insert(0, 1) #The N- and C- cap residues are part of the 1st and last residue index. If no N- and C- cap residues for the protein, comment the line below using "#" del atomsperAminoAcid[1]; del atomsperAminoAcid[-2] # Create Universe uni = mdana.Universe(ref_structure,traj) n,t = list(enumerate(uni.trajectory))[0] box = t.dimensions[:6] atom_Groups = [[] for x in range(chains)] m_start=0 for m in range(0,chains): m_end = atoms * (m+1) atom_Groups[m].extend([uni.select_atoms('bynum '+ str(m_start) + ':' + str(m_end))]) m_start = m_end + 1 fileout1 = open('oligomer-groups.dat','w') fileout2 = open('oligomer-states.dat','w') for tt in uni.trajectory[1:]: fileout1.write (str(tt.frame) + '\t') fileout2.write (str(tt.frame) + '\t') mySet = set([]) graph = [] atom_Groups = [[] for x in range(chains)] m_start=0 for m in range(0,chains): m_end = atoms * (m+1) atom_Groups[m].extend([uni.select_atoms('bynum '+ str(m_start) + ':' + str(m_end))]) m_start = m_end + 1 count = 0 for i in range(chains-1): for j in range(i+1,chains): dist = distance_array(atom_Groups[i][0].positions,atom_Groups[j][0].positions,box).min() if(dist <= Min_Distance): my_tuple = (i,j) mySet.add(my_tuple) graph = nx.from_edgelist(mySet) for i in range(chains): if i not in list(graph): fileout1.write ('['+ str(i)+']' + '\t') fileout2.write ('1' + '\t') else: fileout1.write (str(list(nx.node_connected_component(graph, i))) + '\t') fileout2.write (str(len(list(nx.node_connected_component(graph, i)))) + '\t') fileout1.write ('\n') fileout2.write ('\n') fileout1.close() fileout2.close() # Get oligomerization data OligoStates = [[0 for z in range(chains)] for x in range(frames+1)] file = open("oligomer-groups.dat",'r') line = file.readline() j = 0 while line: temp = line.split('\t') for k in range(chains): OligoStates[j][k] = temp[k + 1][1:-1].split(',') j += 1 line = file.readline() file.close # Create contact matrix ContactMap = [[0 for x in range(AminoAcids)] for y in range(AminoAcids)] # Create atom groups for each amino acid of each monomer AtomGroups = [[] for x in range(chains)] for m in range(0,chains): for aa in range(0,AminoAcids): AtomGroups[m].extend([uni.select_atoms('bynum '+str(atoms*m + atomsperAminoAcid[aa])+':'+str(atoms*m + atomsperAminoAcid[aa + 1] - 1 ))]) count = 0 for n,t in enumerate(uni.trajectory[1:]): on = 0 Groups = [] for i in OligoStates[n]: if len(i) > 1: on = 1 Groups.extend([i]) Set = set(tuple(x) for x in Groups) Groups = [ list(x) for x in Set ] if on == 1: # Calculate dimension of the box to considered PBC box = t.dimensions[:6] for g in Groups: # Calculate contacts for n1,i in enumerate(g): for j in g[(n1 + 1)::]: count += 1 for n2,atoms1 in enumerate(AtomGroups[int(float(i))]): for n3,atoms2 in enumerate(AtomGroups[int(float(j))]): if ((distance_array(atoms1.positions,atoms2.positions,box).min()) <= Min_Distance): ContactMap[n2][n3] +=1 ContactMap[n3][n2] +=1 #print(count) Norm_ContactMap = np.true_divide(ContactMap,float(count)*2.0) # Save contact map in a file fileout = open ('contact-map.dat','w') for i in Norm_ContactMap: for j in i: fileout.write (str(j) + '\t') fileout.write ('\n') fileout.close() #Highest Oligomer size in each frame states=open('oligomer-states.dat', 'r') ter=states.readlines()[0:frames+1] result=[] for freq in (ter): result.append([int(hist) for hist in freq.strip().split('\t')[1:]]) fileout3 = open ('oligo-highest-size.dat', 'w') for oli_count in range(len(ter)): fileout3.write("{} {} {}\n".format(oli_count, '\t', np.max(result[oli_count]))) fileout3.close() # Block Average size_data = np.loadtxt('oligomer-states.dat') window = 25 #over how many frames the running average is to be calculated weights = np.repeat(1.0,window)/window size_data_m = np.convolve(size_data[:,1],weights,'valid') fileout4 = open('oligo-block-average.dat', 'w') for t,b in enumerate(size_data_m): fileout4.write("{} {} {}\n".format(t, '\t', b)) fileout4.close()

6 Appendix 3: Tcl Scripts for the Transition Network Analysis

# Procedure - protein Compactness proc NPMI {chs fr} { set sel [atomselect top "chain $chs" frame $fr] set eig [lsort -increasing -real [lindex [measure inertia $sel eigenvals] 2]] set shape [expr round([lindex $eig 0]/[lindex $eig 2]*10)] $sel delete return $shape } # Procedure - Amino acids in beta-sheet conformation proc beta {chs fr} { set beta [atomselect top "chain $chs and name CA and (betasheet or sheet or beta_sheet or extended_beta or bridge_beta)" frame $fr] set nb [expr [$beta num]] $beta delete return $nb } # Main code # Read input files set input1 [lindex $argv 0] set input2 [lindex $argv 1] # Load trajectory mol new $input1 animate delete beg 0 end 0 animate read xtc $input2 waitfor all # Select peptide length set pep [lindex $argv 2] set pepno [lindex $argv 3] # Define chain id for each protein set ch 0 set nAA [expr $pep*$pepno] set TS {} [atomselect top all] set chain 0 # Assign chain IDs for {set i 0} {$i<$nAA} {incr i $pep} { [atomselect top "residue $i to [expr $i+$pep-1]"] set chain $ch incr ch } # Define cutoff distance set cutoff 4 set d $cutoff set oligomers {} # Calculate the number of frames set nf [molinfo top get numframes] puts "nf = $nf" # Open the transition matrix and states attributes files for writing set fil1 [open "Transition-matrix.dat" w] set fil2 [open "State-Attributes.csv" w] set states {} set prevolig {} set S_unique {} # Cycle through frames for {set j 0} {$j<$nf} {incr j} { set cnt 0 # Go to a specific frame puts "frame $j" animate goto $j display update mol ssrecalc top # Initialize some variables set oligomer {} set oldolig {} set olig {} set transition {} # Cycle through monomers for {set i 0} {$i<$pepno} {incr i} { # Check if the current peptide is already part of the current oligomer if {[lsearch $oldolig $i] == -1} { set cnt {} lappend cnt $i # Identify neighboring chains within distance cutoff set res1 [atomselect top "(within $d of chain $i) and not chain $i"] set length [llength [$res1 get chain]] set neighbor [lindex [lsort -unique [$res1 get chain]] 0] $res1 delete while {$length != 0} { lappend cnt $neighbor set res [atomselect top "(within $d of chain $cnt) and not chain $cnt"] set length [llength [$res get chain]] set neighbor [lindex [lsort -unique [$res get chain]] 0] $res delete } set oldolig [concat $oldolig $cnt] # Update the list of peptides that belong to the same oligomer lappend olig $cnt # Update the list with the oligomer sizes lappend oligomer [llength $cnt] } } # Identify transitions if { $cnt >= 0 } { set n 1 foreach o $olig { # Search for oligomer within the previous oligomer list set so [lsearch $prevolig $o] if {$so>=0} { lappend transition [list [lindex $prevolig $so] $o] } elseif {[llength $o] == 1} { # Search for monomers within previous oligomers (oligomer broken into monomers) foreach p $prevolig { set so1 [lsearch $p $o] if {$so1>=0} { lappend transition [list $p $o] } } } elseif {[llength $o] > 1} { # If the oligomer is larger than 1 look for individual peptides in the previous oligomer list (oligomer formed) foreach ol $o { set so2 [lsearch $prevolig $ol] if {$so2>=0} { lappend transition [list [lindex $prevolig $so2] $o] } else { # Search for individual peptides within previous oligomers foreach pl $prevolig { set so3 [lsearch $pl $ol] if {$so3>=0} { lappend transition [list $pl $o] } } } } incr n } } set prevolig $olig set oldoligomer $oligomer set a2 {} # For each transition identify the aggregation states foreach t1 [lsort -unique $transition] { set a1 {} set b1 {} set S {} foreach t2 $t1 { set frame [expr $j-$f] animate goto $frame mol ssrecalc top lappend a1 [llength $t2] set ss {} set O [llength $t2] # set be 0 #Monomer, call procedures if {$O == 1 } { set Sh [NPMI $t2 $frame] set be [beta $t2 $frame] } else { #Oligomer, call procedures set Sh [NPMI $t2 $frame] set be [beta $t2 $frame] } # Assign a state consisting of three order parameters set states [join $O|[expr round($be/$O)]|$Sh ""] lappend S $states } lappend a2 $a1 set l [lindex $a1 0] set k [lindex $a1 1] if {$cnt!=0} { lappend TS $S } } } } # Write transition matrix and states attributes to files if {1} { set S_unique [lsort -unique [join $TS]] set bbins [llength $S_unique] for {set i 0} {$i < $bbins} {incr i} { for {set j 0} {$j < $bbins} {incr j} { set b($i,$j) 0 } } foreach trans $TS { set i [lsearch $S_unique [lindex $trans 0]] set j [lsearch $S_unique [lindex $trans 1]] set b($i,$j) [expr $b($i,$j)+1] } set row2 {} for {set i 0} {$i < $bbins} {incr i} { set row2 {} for {set j 0} {$j < $bbins} {incr j} { lappend row2 $b($i,$j) } puts $fil1 $row2 puts $row2 } # Create attributes file set all_states {} set count 0 foreach v $TS { if {$count < $pepno} { lappend all_states [lindex $v 0] lappend all_states [lindex $v 1] } else { lappend all_states [lindex $v 1] } incr count } set id 1 puts $fil2 "id state oligomer beta-sheet compactness population" foreach val $S_unique { puts $fil2 "$id $val [lindex [split $val ""] 0] [lindex [split $val ""] 2] [lindex [split $val ""] 4] [llength [lsearch -all $all_states $val]]" incr id } } close $fil1 close $fil2 exit convert2csv.tcl: for converting the transition matrix to csv format #!/usr/bin/tclsh proc splitby { string spl_str } { set lst [split $string $spl_str] for { set cnt 0 } { $cnt < [llength $lst] } { incr cnt } { if { [lindex $lst $cnt] == "" } { set lst [lreplace $lst $cnt $cnt] incr cnt -1 } } return $lst } set input [lindex $argv 0] set output1 [lindex $argv 1] set fil1 [open $input r] set fil2 [open $output1 w] array unset a set i 1 set firstR {} # Read input matrix into variable "a" while {[gets $fil1 line1] >=0} { set firstR [concat $firstR ";$i"] set row [splitby $line1 " "] set j 1 foreach r $row { set a($i,$j) $r incr j } incr i } set n [expr $i-1] # Write output matrix to file puts $fil2 $firstR for {set i 1} {$i <= $n} {incr i} { set var $i for {set j 1} {$j <= $n} {incr j} { set var [concat $var ";$a($i,$j)"] } puts $fil2 $var } close $fil1 close $fil2 quit