field [144]. The protein structure was first minimized with steepest
descent followed by conjugate gradient for 500 steps total. The
system was solvated with 6661 solvent molecules with a NaCl
concentration of 0.2 M in water using the TIP3P water model.
The final system was minimized once more for 2000 steps to
remove bad contacts. Consecutive position-restrained 1 ns canonical ensemble and 1 ns isothermal–isobaric ensemble simulations
were performed with a heavy atom restraint force constant of
500 kcal/(mol A ˚ 2 ) with the Berendsen thermostat and barostat
[145] at 300 K and 1 bar, respectively. The system was then
simulated under isothermal–isobaric ensemble for 700 ps converging the system density, followed by a 5 ns production simulation.
System snapshots were taken at 100 ps intervals along the 5 ns
trajectory and were further simulated using the microcanonical
ensemble (NVE) for 150 ps using 0.5 fs time steps, outputting
the velocity and position trajectory files every 10 fs.
The trajectories of the MD simulations were then analyzed. To
compute the inter-residue flux for each NVE simulation we used
the CURP (CURrent calculations for Proteins) program developed
by Yamato and coworkers [13]. The flux between residues is multiplied by k B T, where T was 300 K in the MD simulations, and the
results are reported in units of (kJ/mol)
2 ps
À1 . The polar contact
dynamics was calculated over the same trajectory. A polar contact
search was carried out on each simulation looking for all donor
nitrogen or oxygen, donor proton and oxygen acceptor triplets and
for all triplets where the donor proton and oxygen acceptor distances are 0.28 nm for at least 99% of the NVE trajectory. For the
pairs that met the distance criteria we calculated 1= δr
2
ij
D
E
. For each
trajectory the polar contact dynamics calculated was paired with the
corresponding energy flux. We then defined hydrogen bonds as
having a donor-proton-acceptor angle greater than or equal to
150
over the duration of the simulation.
For each pair found in any one of the simulations that satisfied
the hydrogen bond criterion, we plot the value of the flux against
1= δr
2
ij
D
E
in Fig. 6. As in the HP36 work [23], amino acid pairs
within 4 in sequence space are not included, since those pairs are
within helices along which energy transport is anyway relatively fast,
the rate constants for which follow a different relation appropriate
for energy transfer along the backbone [23]. The data are separated
into two groups, one for hydrogen bond pairs that lie 5–9 amino
acids away and another group where they are further away in
sequence space. The first group consists only of pairs that are
5 and 6 away in sequence space, and 3 of the four pairs identified
involve backbone–backbone hydrogen bonds. They are shown in
Fig. 7, and are seen to be at turns just beyond one of the helices.
The points corresponding to those hydrogen bonds appear to be
well fit by a line in Fig. 6 given by flux ¼ 76.4 þ 0.0273/ δr
2
ij
D
E
,
50
Korey M. Reid and David M. Leitner
descent followed by conjugate gradient for 500 steps total. The
system was solvated with 6661 solvent molecules with a NaCl
concentration of 0.2 M in water using the TIP3P water model.
The final system was minimized once more for 2000 steps to
remove bad contacts. Consecutive position-restrained 1 ns canonical ensemble and 1 ns isothermal–isobaric ensemble simulations
were performed with a heavy atom restraint force constant of
500 kcal/(mol A ˚ 2 ) with the Berendsen thermostat and barostat
[145] at 300 K and 1 bar, respectively. The system was then
simulated under isothermal–isobaric ensemble for 700 ps converging the system density, followed by a 5 ns production simulation.
System snapshots were taken at 100 ps intervals along the 5 ns
trajectory and were further simulated using the microcanonical
ensemble (NVE) for 150 ps using 0.5 fs time steps, outputting
the velocity and position trajectory files every 10 fs.
The trajectories of the MD simulations were then analyzed. To
compute the inter-residue flux for each NVE simulation we used
the CURP (CURrent calculations for Proteins) program developed
by Yamato and coworkers [13]. The flux between residues is multiplied by k B T, where T was 300 K in the MD simulations, and the
results are reported in units of (kJ/mol)
2 ps
À1 . The polar contact
dynamics was calculated over the same trajectory. A polar contact
search was carried out on each simulation looking for all donor
nitrogen or oxygen, donor proton and oxygen acceptor triplets and
for all triplets where the donor proton and oxygen acceptor distances are 0.28 nm for at least 99% of the NVE trajectory. For the
pairs that met the distance criteria we calculated 1= δr
2
ij
D
E
. For each
trajectory the polar contact dynamics calculated was paired with the
corresponding energy flux. We then defined hydrogen bonds as
having a donor-proton-acceptor angle greater than or equal to
150
over the duration of the simulation.
For each pair found in any one of the simulations that satisfied
the hydrogen bond criterion, we plot the value of the flux against
1= δr
2
ij
D
E
in Fig. 6. As in the HP36 work [23], amino acid pairs
within 4 in sequence space are not included, since those pairs are
within helices along which energy transport is anyway relatively fast,
the rate constants for which follow a different relation appropriate
for energy transfer along the backbone [23]. The data are separated
into two groups, one for hydrogen bond pairs that lie 5–9 amino
acids away and another group where they are further away in
sequence space. The first group consists only of pairs that are
5 and 6 away in sequence space, and 3 of the four pairs identified
involve backbone–backbone hydrogen bonds. They are shown in
Fig. 7, and are seen to be at turns just beyond one of the helices.
The points corresponding to those hydrogen bonds appear to be
well fit by a line in Fig. 6 given by flux ¼ 76.4 þ 0.0273/ δr
2
ij
D
E
,
50
Korey M. Reid and David M. Leitner