Thermostat

MD simulations were performed using the models generated by homology as the initial coordinates. Set up of the systems was performed with GROMACS 4.5 [85 (link)]. First, each system comprising in about 200,000 atoms was energetically relaxed in vacuum using the steepest descend combined with the conjugate gradient algorithm, then solvated in a TIP3P water box and neutralized by adding counter-ions (Cl^-), resulting in 200 000 atoms systems. Each solvated system was relaxed by 30 000 steps of minimization using the conjugate gradient algorithm. In the first 10 000 steps of minimization, constraints were applied to the protein heavy (non-hydrogen) atoms. In the following 10 000 steps, each system was minimized with constraints applied on Cα atoms only while in the last 10 000 steps, no constraints were applied. After relaxation, each system was linearly heated from 10 to 310 K with constraints on Cα. An unconstrained MD simulation was then performed at constant volume (NVT) using a Berendsen thermostat [86 (link)] for 100 ps and for further 100 ps at constant pressure (1 bar) using the Parrinello-Rahman algorithm [87 (link)]. A last MD simulation at 310 K and 1 bar (NPT ensemble) were carried out for 5 ns to achieve the properly equilibrated models.

The β₂-AR-agonist-NDP complex, with peptide coordinates obtained from the docking output, was inserted into a 70 × 70 Å 2-oleoyl-1-palmitoyl-sn-glyecro-3-phosphocholine (POPC) lipid bilayer. A 10 Å water layer was added from the positive side of the z-axis and 60 Å from the negative side. Bad-contact water molecules were removed from the lipid membrane bilayer using the appropriate tcl script. Additionally, the system was neutralized with 0.15 M NaCl, resulting in a 61,183 (~60,000) atom ensemble. The system was subject to a 10,000 step energy minimization, 250 NVE ps equilibration, and 100 ns NPT MD production. Pressure and temperature were set to 1 bar and 310 K, respectively, using a Berendsen thermostat, and the applied integration step was 1 fs. In all simulations, periodic boundary conditions with particle-mesh Ewald calculations were implemented. The cut-off was set to 12 Å. A CHARMM22^{22 (link),23 (link)} force field was used for protein and lipids, and CGenFF^{43 (link),44 (link)} was used for ligands.

The lipid bilayer was composed of POPC, POPE and tetraoleoyl cardiolipin in a ratio of about 9:7:4, mimicking the average lipid composition of inner mitochondrial membranes from guinea pig liver, rat liver and pig heart [65 (link)]. An initial bilayer was constructed from 125 lipid molecules (56 POPC, 44 POPE, and 25 cardiolipins) by first placing lipids randomly in a simulation box, with six times as many coarse-grained water beads. Lipid bilayers then self-assembled during a simulation of 60 ns, with pressure coupled isotropically using a Berendsen barostat (coupling constant of 3 ps), and temperature with a Berendsen thermostat (coupling constant of 0.3 ps). The self-assembled bilayer was then simulated for a further 15 ns. A larger bilayer of 250 lipids was assembled by pasting together two of the smaller self-assembled bilayers. A second bilayer model was generated by simulating the first bilayer for a further 300 ns.

The GROMACS software package [36 (link)] was used to perform MD simulations. The GROMOS force-field 43a2 [37 ] was used to describe the peptide and peptide–solvent interactions. The force-fields were parametrized for use with a group-based twin range cut-off scheme (using cutoffs of 1.0/1.4 nm and a pair-list update frequency of once per 10 steps), including the particle-mesh Ewald (PME) method. The water was modeled using the SPC model [38 (link)]. A time step of 2 fs was used. Bond lengths were constrained using the LINCS algorithm [39 (link)]. The simulations were performed in the NPT ensemble using periodic boundary conditions. The temperature was weakly coupled (coupling time 0.1 ps) to T = 298 K using the Berendsen thermostat [40 (link)]. The pressure was also weakly coupled (coupling time of 1.0 ps and compressibility of 4.5 × 10⁻⁵), using an isotropic coupling scheme at 1 bar. The DSSP protocol and the corresponding computer program [41 (link)] were employed to analyze the time evolution of the local secondary structure.

For each equilibrated system (models I–V) two independent MD simulations were run with different initial velocities using the PMEMD module of AMBER 10. The temperature was kept at 310 K (Berendsen thermostat), and pressure at 1 bar (Langevin piston coupling algorithm). The SHAKE algorithm was used to freeze covalent bonds involving hydrogen atoms, allowing for an integration time step of 2.0 fs. Long-range electrostatic interactions were treated by the Particle Mesh Ewald method [54] . Each simulation run was continued until a total time of 50 ns to 70 ns, depending on the estimated convergence of the first 50 ns of simulation. Coordinates files were recorded every 1 ps.

All non-bonded interactions were shifted to zero when the distance between beads was 12 Å. A timestep of 20 fs was used and the temperature was held at 310 K using a Berendsen thermostat with a time constant of 1 ps. The pressure was maintained at 1 bar by a Berendsen barostat, applied semi-isotropically, with a time constant of 1 ps and a compressibility of 3 × 10^–4 bar.²² Coordinates of all molecules were written to disc every 2 ns.