*3.9. Molecular Modelling*

Molecular modelling calculations were performed on E4 Server Twin 2× Dual Xeon-5520, equipped with two nodes. Each node: 2× Intel® Xeon® QuadCore E5520-2.26Ghz, 36 GB RAM. The molecular modelling graphics were carried out on a personal computer equipped with Intel(R) Core(TM) i7-4790 processor and SGI Octane 2 workstations.

*Conformational property analysis.* The apparent pKa and logD values (pH 7.4 and 7.2) of compounds were calculated by using the ACD/Percepta software (ACD/Percepta, Advanced Chemistry Development, Inc., Toronto, ON, Canada, 2017, http://www.acdlabs.com).

The compounds were considered neutral in all calculations performed as a consequence of the estimation of percentage of neutral/ionized forms computed at pH 7.4 (blood pH value) and pH 7.2 (cytoplasm pH value) using the Henderson–Hasselbalch equation. The compounds were built using the Small Molecule tool of Discovery Studio 2017 (Dassault Systèmes BIOVIA, San Diego). Then, the built compounds were subjected to molecular mechanic (MM) energy minimization (ε = 80 × *r*) until the maximum RMS derivative was less than 0.001 kcal/Å, using Conjugate Gradient [58] as minimization algorithm. Atomic potentials and charges were assigned using the CFF forcefield [48]. The conformers obtained for each compound were used as starting structure for the subsequent systematic conformational analysis (Search Small Molecule Conformations; Discovery Studio 2017). The conformational space of the compounds was sampled by systematically varying the rotatable bonds sp3–sp3 and sp3–sp2 with an increment of 60◦. The RMSD cutoff for structure selection was set to 0.01 Å. Finally, to ensure a wide variance of the input structures to be successively fully minimized, an energy threshold value of 10<sup>6</sup> kcal/mol was used as selection criteria. The generated structures were then subjected to MM energy minimization (CFF forcefield; ε = 80 × *r*) until the maximum RMS derivative was less than 0.001 kcal/Å, using Conjugate Gradient as minimization algorithm. Finally, the resulting conformers were ranked by their potential energy values (i.e., ΔE from the global energy minimum (GM)).

All MM conformers within 5 kcal/mol from GM has been then subjected to DFT calculations. The calculations were carried out using the Gaussian 09 package [59]. All structures were fully optimized at the B3LYP/6-31+G (d,p) level using the conductor-like polarizable continuum model (C-PCM) [49,60,61]. In order to characterize every structure as minimum and to calculate the Gibbs free energy a vibrational analysis was carried out at the same level of theory using the keyword: freq. The RMS force criterion was set to 3 × 10−<sup>4</sup> a.u. Partial charges, molecular orbitals and spin density have been calculated using the natural bond orbital (NBO) method [62]. The resulting conformers were ranked by their potential energy values (i.e., ΔE from the global energy minimum (GM)) and classified by their dihedral angle values).

*Calculation of redox properties.* An appropriate way for calculating the redox potential is by using a thermodynamic cycle. Accordingly, the redox potentials for the compounds was calculated by using the Born–Haber cycle (Scheme 2), which links the Gibbs free energy change of the one-electron transfer reaction in the gas phase with that of the same reaction in aqueous solution [63,64].

$$\begin{array}{ccc} \mathsf{X}\_{\mathsf{(g)}} + \mathsf{e}^{\cdot} \to \mathsf{X}\_{\mathsf{(g)}}^{\bullet} & \mathsf{\Delta G}^{\circ} \text{ (gas)} \\\\ \mathsf{\Delta G}^{\circ} \times \mathsf{X}\_{\mathsf{(z\lhd\mathsf{v})}} \Big| & \mathsf{\Delta G}^{\circ} \mathsf{X}\_{\mathsf{(z\lhd\mathsf{v})}} \\\\ \mathsf{X}\_{\mathsf{(z\lhd\mathsf{j})}} + \mathsf{e}^{\cdot} \to \mathsf{X}\_{\mathsf{(z\lhd\mathsf{j})}}^{\bullet} & \mathsf{\Delta G}^{\circ} \text{ (red)} \end{array}$$

**Scheme 2.** Born–Haber cycle for a generic one-electron transfer reaction in vacuo and in aqueous solution.

Accordingly, the Born–Haber thermodynamic cycle allows to include the desolvation/solvation effects by calculating the Gibbs free energy values of the reaction in gas phase and in solution. According to this approach, the standard Gibbs free energy of the electron transfer ( Δ *<sup>G</sup>*0*red*,*aq*), was calculated taking into account the free energy change in the gas phase ( Δ *<sup>G</sup>*0*gas*) and the solvation free energies of the oxidized ( Δ *<sup>G</sup>*0*solv*(*X*)) and reduced species Δ *<sup>G</sup>*0*solv*(*X*•−)) as reported in the following equation:

$$
\Delta G^{0}\_{\text{red,aq}} = \Delta G^{0}\_{\text{gas}} + \left(\Delta G^{0}\_{\text{solv}}(X^{\bullet -}) - \Delta G^{0}\_{\text{solv}}(X)\right) \tag{1}
$$

Accordingly, we, firstly, calculated the Gibbs free energy in gas phase by using the following equation:

$$
\Delta G^{0}\_{\text{gas}} = G^{0}\_{\text{gas}}(X^{\bullet-}) - G^{0}\_{\text{gas}}(X) \tag{2}
$$

Then, the Gibbs free energy of solvation of both the species X and X•− were calculated according to the following equations:

$$
\Delta \mathcal{G}^{0}\_{solv}(X) = \mathcal{G}^{0}\_{aq}(X) - \mathcal{G}^{0}\_{\text{gas}}(X) \tag{3}
$$

$$
\Delta G\_{\text{solv}}^{0}(X^{-}) = G\_{a\eta}^{0}(X^{\bullet-}) - G\_{\text{gas}}^{0}(X^{\bullet-}) \tag{4}
$$

The values obtained were used to calculate the standard Gibbs free energy of the overall reaction (Δ *<sup>G</sup>*0*red*,*aq*; kcal·mol−1) according to Equation (1). The X and X•− species correspond in our case to Q and Q•−, respectively in the first electron transfer, and to QH• and QH− in the second electron transfer (see Scheme S1).

Finally, the Δ *<sup>G</sup>*0*red*,*aq*, were used to calculate the standard reduction potential by the Nernst equation:

$$E^0 = -\frac{\Delta G^0\_{red,aq}}{nF} \tag{5}$$

where *n* is the number of electrons involved in the process (i.e., 1) and *F* the Faraday constant (23.06 kcal mol−<sup>1</sup> <sup>V</sup>−1).

The reduction potential given by the Equation (5) is an absolute value of the parameter and has to be referred to the standard hydrogen electrode (SHE; redox potential of 4.43 V), to obtain the standard redox potential. Accordingly, the value of 4.43 V was subtracted from the values obtained by Equation (5).

On these bases, in order to calculate the redox potential, starting from the structure of the DFT Q GM conformers (i.e., the starting quinone species), the redox states Q•−, QH•, QH<sup>−</sup>, QH2 were generated. Following the reduction pathway of quinones showed in Scheme S1, each species was generated starting from the DFT optimized species of the previous step and submitted to full DFT optimization both in aqueous solution (using the C-PCM method) and in vacuo using the same parameters above described.

The standard Gibbs free energy for the protonation of the reduced species were calculated by determining Δ G◦ of the reaction Y + H<sup>+</sup> → YH that is Δ G◦ = G◦(YH) − G◦(Y), where Y is Q•− or QH− DFT conformer. The energy of the frontier molecular orbitals (HOMO, LUMO, and SOMO) were considered for the species Q, Q•−, QH• and QH− fully optimized by DFT calculations performed in the

aqueous solution (using the C-PCM method). According to the Koopmans' theorem [65], the ionization potential (IP) is given by the negative of the value of the energy of the HOMO.
