*4.3. Molecular Dynamics Simulations*

The CHARMM 36 force field [34] was used for the simulations summarized in Table 1 using the GROMACS software package version 5.1.3 [35]. The YqhD dimer bound to Zn2+ and NADP/H cofactors was centered in a cubic periodic box (~11 nm<sup>3</sup> ) and set to have a distance larger than 1 nm from any side of the box. Solvent molecules having any atom within 0.15 nm from the protein were removed. TIP3P model [36,37] was used for explicit water solvent. Sodium counter ions were added by replacing the solvent molecules at the sites of most negative electrostatic potential to provide the box with a total charge of zero. The protonation state of residues was assumed to be the same as that of the isolated amino acids in solution at pH 7. The LINCS [38] algorithm was used to constrain all bond lengths and the SETTLE [39] algorithm was used for the water molecules. Electrostatic interactions were calculated using the Particle Mesh Ewald method [40]. For the calculation

of long-range interactions, a grid spacing of 0.12 nm combined with a fourth-order B-spline interpolation was used to compute the potential and forces between grid points. A nonbonded pair-list cutoff of 1.4 nm was used and updated at every five time-steps. V-rescale thermostat [41] was used to keep the temperature at 300 K through a weak coupling of the system to an external thermal bath with a relaxation time constant τ = 0.1 ps. The pressure of the system was kept at 1 bar using Berendsen's barostat [42] with a time constant of 1 ps. A time step of 2 fs was used to integrate the equations of motion.


**Table 1.** Simulation summary for the YqhD homodimer in aqueous solution.

First, the simulated systems were energy minimized, using the steepest descent algorithm, for at least 5000 steps to remove clashes between atoms that were too close. After energy minimization, all atoms were given an initial velocity obtained from a Maxwell-Boltzmann velocity distribution at 300 K to start the MD simulations. The system was initially equilibrated by 30 ps with position restraints on the heavy atoms of the dimer to allow relaxation of the solvent molecules. After the equilibration procedure, position restraints were removed, and the system was gradually heated from 50 K to 300 K during 200 ps of simulation. The equilibrated structure was used to perform MD simulations in the NPT ensemble for 20 ns at 300 K. The final conformation of the 20 ns isothermalisobaric NPT ensemble was used to set up five sets of 200 ns independent simulations in the canonical NVT ensemble initiating with five different velocities.

The starting coordinates were adopted from the crystal structure which has modified cofactor (NADPH(OH)2); hence, the equilibrated conformations of YqhD bound to NADP and NADPH cofactors were obtained after a 20 ns NPT simulation, which was used as the starting point for the final production run of 200 ns long NVT simulations. The starting crystallographic coordinates of the YqhD homodimer were used as the reference structure for the analysis of the trajectories to characterize the conformational changes induced by a change in the oxidation state of NADP/H cofactor. The analysis was focused on three sets/sub-structures within the simulation data, considering dimer, monomer, and domain in a 1 µs aggregated trajectory for each cofactor oxidation state (NADP and NADPH). Structural analysis includes the root mean square deviation (RMSD), root mean square fluctuation (RMSF), radius of gyration (Rg), and cofactor binding interactions (Zn+2 and NADP/H), comparing these with respect to the crystal structure.
