**2.3 Molecular coupling**

*Homology Molecular Modeling - Perspectives and Applications*

previous studies.

NPT dynamics at 300 K and 1 bar.

**2.1 Gaussian accelerated molecular dynamics (GaMD)**

each one representing a different energy and conformational state.

The selected nanomaterials were Au nanoclusters (AuNCs) modified with two types of linkers: 4-aminothiophenol (ATP) and mercaptobenzoic acid (MBA); each cluster had 2.1 nm diameter and a nucleus composed of 96 Au atoms (Au314 (SR) 96). The nanostructures were modeled as reported elsewhere [54]. Six gold

**2.2 Building of nanomaterial structures**

**2. General techniques: Horseradish modeling and system building**

The X-ray crystal structure from horseradish peroxidase (HRP) was collected from Protein Data Bank (www.rcsb.org/, PDB code 1H5A) [45], and reconstructed using i-tasser metaserver [46]. The best model was selected based on lowest RMSD and higher C-value, to achieve the best initial coordinates. To improve the model system, all missing atoms and sidechain positions were modeled against available plant peroxidases on the PDB. A penta-coordinated resting state heme group was crystallized as Fe3+, and axial coordination with a histidine groups was retained. The Protein Preparation Wizard routine (Schrödinger Maestro v2019-4, New York, 2019), was applied as a preparation step to correctly assign missing hydrogen atoms and protonation states of ionizable residues. Finally, PROPKA2.0 sub-routine was applied to whole model system [47] and pH 7 was fixed accordingly to our

Enzyme structural parameters were described using the ff14SB force field [48]. The full systems comprised ~60,000 atoms in a cubic box of 15 Å length including TIP3P explicit water and Na + ions to ensure overall charge neutrality. Non-bonded interactions were calculated within a 12 Å cutoff, and long-range electrostatics were treated using the Particle-Mesh Ewald method [49]. SHAKE algorithm was enabled to constrain all bonds involving hydrogen during simulations. Conventional MD protocol for each system replica comprised: 1) 5000 steepest descent minimization steps followed by 10,000 conjugate gradient minimization steps; 2) 500 ps of progressive NVT heating from 0 to 300 K; 3) 5 ns of NVT equilibration, and 4) 10 ns of

Conventional MD simulations were carried out to initially relax full system coordinates, using AMBER18 software package [50]. Then, the GaMD module implemented in AMBER v. 18 was applied to perform extra 50 ns of short classic MD simulation to collect the statistics for calculating GaMD acceleration parameters. Extra 50 ns of short equilibration was applied after adding the boost potential, and finally three independent 500 ns GaMD production simulations with randomized initial atomic velocities. This method applies Gaussian functions (boost potentials) on the total potential energy term, which disturb the potential energy surface allowing the exploration of different energy states. In addition, functions that directly disturb the dihedral angles of the amide bonds were applied, promoting conformational changes. All GaMD simulations were performed with a dual-boost level by setting the reference energy to the lower bound. One boost potential is applied to the dihedral energetic term and the other to the potential energy term. The average and SD of the system potential energies were calculated every 500,000 steps (1 ns) for all simulation systems. The upper limit of the boost potential SD, σ0, was set to 6.0 kcal·mol−1 for both the dihedral and the total potential energetic terms [51–53]. The system temperature was ~298 K and 1 atm pressure, with integration steps of 2 fs. Three different system coordinates were extracted from GaMD trajectories,

**106**

Adaptive Steered Molecular Dynamics (ASMD) simulations [55, 56] started from the earlier selected conformational states, derived from the GaMD analysis. Due to the chemical complexity of the enzyme surface, the simulations were driven by the electrostatic surface (PME) [49] for each conformational state and the 6 AuNCs layer. Hence, the scans were performed at speeds of 10 Å / ns, which comprised a 10-step profile of 1 Å between the Fe(III) of the active site and the ATP amino group or MBA group carboxylic of the 6 AuNCs. Three independent replicas were performed through gradual scans, involving measurements of the free energy of coupling between HRP and 6 AuNCs. The free energy was measured by shortening the distance between the central Fe(III) atom and the amino groups of the linkers, i.e., ATP amino groups and central Fe(III) atom of HRP heme (NH3 + ➔ Fe3+).
